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1.  Introduction 


It  is  well  recognized  that  laminated  composite  plates  exhibit  extraordinary  strength 
per  weight  ratio  characteristics.  Consequently,  their  use  in  aerospace  structures  is  the  way 
of  the  future.  However,  in  order  to  be  able  to  predict  their  failing  characteristics 
particularly  in  the  neighborhood  of  free  surfaces  such  as  straight  edges,  holes  edges,  cracks 
etc.,  it  is  necessary  to  know  the  local  stress  behavior  particularly  from  a  3D  point  of  view. 
Thus,  if  rational  designs  using  fiber  reinforced-resin  matrix  composite  laminates  are  to  be 
made,  their  performance  under  static,  dynamic  and  fatigue  loads  need  to  be  predictable. 
The  first  step  towards  this  goal  is  the  realization  that  the  ultimate  failure,  as  well  as  many 
other  aspects  of  the  composite  behavior,  is  the  result  of  the  growth  and  accumulation  of 
micro-damage  to  the  fibers,  matrix  and  their  interfaces.  Thus,  it  appears  that  any  generally 
successful  model  of  performance  and  failure  must  incorporate  the  effects  of  this  local 
damage  in  some  way.  This  certainly  represents  a  challenge. 

In  this  investigation,  a  systematic  ,  3D,  micromechanics  approach  is  used  where 
the  fibers  are  assumed  to  be  cylindrical  inclusions  which  are  embedded  into  a  matrix  plate. 
A  three  dimensional  analysis  is  used  in  order  to  capture  any  edge  effects  which  may  be 
present.  A  set  of  fundamental  key  problems  is  identified  and  their  respective  solutions  are 
then  used  in  order  to  provide  some  answers  to  the  following  fundamental  questions: 

(i  )  UNIDIRECTIONAL  COMPOSITE  PLATE: 
transverse  strength 
edge  delamination 
modeling  of  fiber  matrix  interface 
longitudinal  strength 

residual  stresses  due  to  thermal  expansion  mismatch 

The  approach  consists  of  first  solving  for  the  3D  stress  field  of  one  fiber  only.  The 
solution  is  then  extended  to  the  case  of  a  periodic  array  of  fibers.  The  3D  results  are  then 
used  to  first  identify  critical  locations  where  failure,  due  to  fracture,  is  most  likely  to  initiate 
and  second  to  derive  fracture  criteria  for  crack  initiation  at  the  local  level.  The  criteria  reveal 
the  dependance  of  the  transverse  strength  on  the  material  properties,  the  local  geometry  and 
the  applied  loads.  Ultimately,  a  general  criterion  will  be  derived  which  incorporates  all  of 
the  above  effects.  In  this  phase  of  the  research  we  are  focusing  on  the  process  of  crack 
initiation  and,  as  a  result,  the  criteria  will  be  valid  only  for  crack  lengths  up  to  30%  or  40% 


of  the  fiber  radius.  This  is  because  the  present  analysis  does  not  account  for  the  interactions 
between  cracks  and  periodic  fibers.  It  may  be  noted,  however,  that  it  is  now  possible  to 
indeed  take  into  account  such  interactions  at  the  local  level,  a  subject  which  will  be 
investigated  in  the  near  future.  Be  that  as  it  may,  the  identification  of  the  critical  areas  for 
possible  failure  and  the  prediction  of  the  critical  stress  for  crack  initiation  are  of  great 
interest  to  the  desingners. 

(ii)  BRIDGING  THE  GAP  BETWEEN  MICRO  AND  MACROMECHANICS 

Our  investigation  shows  that  it  is  now  possible  to  relate  results  based  on 
macTomechanical  considerations  with  results  based  on  micromechanical  considerations  via 
certain  correlation  functions  which  are  obtained  in  part  (i).  As  a  practical  matter ,  once  these 
correlation  functions  have  been  established,  then  results  based  on  macromechanical 
considerations  (eg.  finite  ellements  etc.)  may  then  be  used  to  also  answer  questions  of 
damage  at  the  local  level.  Thus,  one  of  the  important  findings  in  part  (i)  is  the  identification 
of  the  regions  of  applicability  of  the  results  which  are  based  on  macromechanical 
considerations. 

(iii)  LAMINATED  COMPOSITE  PLATE  BASED  ON  MICROMECHANICS 

The  above  work  will  next  be  extended  to  the  more  general  case  of  a  laminated 
composite  plate  of  different  ply  orientations.  The  cylindrical  fibers  in  this  case  have  two 
different  fiber  orientations  ,  thus  modeling  from  a  micromechanics  point  of  view  the 
consept  of  an  interface. 

2a.  Effects  of  Transverse  Load  on  a  Fiber/Matrix  Interface  Failure 

The  low  transverse  tensile  strength  of  unidirectional  laminae  presents  a  major 
problem  in  the  design  of  composite  laminate  structures.  Although  the  fibers  can  be  oriented 
so  that  they  are  parallel  to  the  external  loads,  it  is  almost  impossible  to  avoid  transverse 
stresses  which  may  lead  to  premature  failiu«  of  the  laminate.  An  excellent  example  of  this 
is  in  the  design  of  a  filament  wound  pipe  or  pressure  vessel . 


There  is  no  simple  relation  for  predicting  the  transverse  strength.  Unlike  the 
longitudinal  tensile  strength  which  is  determined  almost  entirely  by  a  single  factor,  i.e.  the 
fiber  strength,  the  transverse  strength  is  governed  by  many  factors  including  the  properties 


of  the  fiber  and  matrix,  the  interface  bond  strength,  the  presence  and  distribution  of  voids, 
and  the  internal  stress  and  strain  distribution  due  to  the  interactions  between  fibers,  voids 
etc. 


Recently,  we  were  able  to  investigate  what  effect  if  any  does  a  transverse  load  have 
on  the  failure  characteristics  of  a  fiber/matrix  interface.  To  establish  this,  we  carried  out 
an  analysis  based  on  3D  analytical  considerations  .  (  For  more  details  see  Part  I).  The 
analysis  shows  the  stresses  and  oee  for  ratios  of  (a/h)  <  0.05,  to  be  almost  constant  in 
the  interior  of  the  laminate  and  that  as  one  approaches  the  free  edges  of  the  laminate,  i.e.  the 
region  where  the  fiber  intersects  the  free  surface,  a  boundary  layer  effect  is  shown  to 
prevail.  A  separate  local  asymptotic  analysis  carried  out  by  the  author  shows  that  the  stress 
field  there  is  indeed  singular.  For  example,  for  a  glass  fiber/  epoxy  matrix  the  stress 
singularity  is  0.25  (  Folias  1989),  while  for  the  case  of  a  carbon  fiber/epoxy  matrix  it  is 
0.31  ( for  more  details  see  Part  II)  .This  neighborhood,  therefore,  is  an  area  with  potential 
for  a  fiber/matrix  interface  failure  due  to  the  relatively  high  stress  levels  which  are  present 
It  may,  furthermore,  be  noted  that  for  ratios  of  (a/h)  >  0.3  the  stress  profiles  along  the 
fiber  length  is  very  much  non  uniform.  The  model,  however,  is  still  applicable  and  it  does 
provide  us  with  accurate  results. 

The  results  may  now  be  used  to  derive  failure  criteria  based  on  micromechanical 
considerations.  For  example,  such  a  criterion  due  to  the  presence  of  an  interface  crack  has 
been  derived  and  the  details  may  be  found  in  PART  III  (see  eq.34  and  fig  12).  Similarly  a 
failure  criterion  due  to  the  presence  of  a  0  -  type  of  crack  has  also  been  derived. 

Finally,  in  the  limit ,  if  one  allows  all  the  fibers  to  debond  completely  one  arrives  at 
a  local  lower  bound  for  the  overall  matrix  strength  of  the  material  system.. 

2b.  Effect  of  Fiber  Coating  on  the  Strength  Characteristics 

The  properties  of  fiber  reinforced  composites  are  very  much  dependendent  upon  the 
interfacial  region  between  the  matrix  and  the  fiber.  This  is  because  the  primary  function  of 
this  interface  region  is  to  transmit  a  portion  of  the  load  from  the  matrix  to  the  reinforcing 
fibers  and  vise  versa.  This  ability  to  transmit  stress  across  the  phase  boundary  must, 
ultimately,  depend  upon  the  mechanical  properties  of  the  matrix,  the  load  bearing  capacity 
of  the  fibers  as  well  as  the  strength  of  the  fiber  matrix  interface.  It  is  natural,  therefore,  to 
seek  the  relationship  between  the  overall  composite  strength  and  the  above  variables. 


Based  on  3D  considerations  the  subject  was  recently  investigated  by  Folias  and  Liu  (  see 
PART  rv )  for  transverse  loads  and  by  Zhong  and  Folias  ( see  PART  VI )  for  axial  loads. 

In  the  case  of  a  glass  fiber/  epoxy  matrix  the  effect  of  a  fiber  coating  is  shown  to 
have  a  minimal  effect  on  the  stress  field.  For  other  type  of  material  systems,  however,  this 
effect  can  be  substantial.  Our  3D  model  does  provide  us  with  accurate  results.  For  more 
details  see  PART  IV.  Finaly,  as  a  practical  matter,  if  one  has  a  functionally  gradient 
property  matrix  then  the  singular  stress  at  the  fiber  edge  may  also  be  eliminated. 

3.  The  Effect  of  Thermal  Stresses  and  Curing  Stresses 

Shrinkage  stresses  during  cure  and  thermal  stresses  due  to  differences  between  the 
thermal  expansion  coefficients  of  the  matrix  and  fiber  may  have  a  major  effect  on  the 
microstresses  within  a  composite  material  and  as  a  result  their  contributions  should  be 
added  to  the  stresses  introduced  by  the  mechanical  loads.  It  may  also  be  noted  that  these 
microstresses  are  often  sufficient  to  produce  microcracking  even  in  the  absence  of  external 
loads. 


As  it  was  noted  previously,  the  stress  field  in  the  neighborhood  where  a  fiber  meets 
a  free  surface  is  singular.  Thus,  all  things  being  equal,  small  cracks  are  more  likely  to 
develop  at  such  neighborhoods  and  subsequently  propagate  into  the  interior  of  the  matrix. 
The  condition,  as  expected,  is  further  aggravated  when  a  temperature  change  takes  place 
because  of  the  difference  in  the  thermal  expansion  coefficients  between  matrix  and  fibers. 
The  problem  is  of  panicular  concern  to  ceramic  type  of  composites. 

The  objective  of  this  investigation  was  to  study  the  3D  stress  field  of  a  periodic 
array  of  fibers  embedded  into  a  matrix  plate  and  subjected  to  a  temperature  change  2/r,  e.g. 
during  the  curing  process.  Particular  emphasis  is  placed  in  the  neighborhood  of  the  free 
surface  as  well  in  the  interior  of  the  fiber  length.  The  analysis  reveals  the  influence  that 
each  individual  ratio,  fiber  radius  to  fiber  length ,  fiber  volume  fraction,  material  properties 
Gf/Gm,  cooling  temperature,  and  thermal  expansion  coefficients  have  on  the  residual 
stresses  field  which  is  induced  within  the  material  system  during  the  process  of 
manufacturing  for  example.  The  model  identifies  the  critical  locations  where  possible 
failures  are  most  likely  to  initiate.  In  general,  the  location  dependents  on  the  material 
properties  as  well  as  the  other  ratios  previously  mentioned.  For  more  details,  see  PART  V. 
Fracture  criteria  for  crack  initiation  have  also  been  derived  and  the  results  will  be  reported 


soon.  It  may  be  noted  that  the  model  is  applicable  to  intermetalic,  ceramic  and  organic 
composites.  The  applications  are  numerous,  e.g.  jet  engine  technology,  automobile 
engines,  thermal  fatigue  for  high  performance  aircraft,  etc. 

4.  Zone  of  Influence  of  a  Break  in  one  Fiber 

Unidirectional  fiber  composites  are  important  structural  elements  of  many  modem 
composite  materials.  Their  tensile  strength  frequently  determines  the  possibility  of  using 
them  in  various  structures.  Due  to  various  technological  reasons,  breakage  of  individual 
fibers  long  before  fracture  of  the  entire  specimen  is  absolutely  unavoidable.  For  example, 
in  glass-reinforced  plastics,  some  individual  fibers  break  at  loads  of  only  one-tenth  of  the 
maximum  load  for  the  material.  Actually,  all  fibers  used  in  composites  usually  have  high 
mean  statistical  strength,  but  an  extremely  low  fracture  toughness,  and  as  a  result,  they  are 
quite  sensitive  to  tiny  cracks,  which  unavoidably  arise  in  the  process  of  manufactme  of  the 
fibers.  Therefore,  a  certain  number  of  these  fibers  always  break  at  extremely  low  tensile 
loads. 


In  view  of  the  above,  the  following  questions  come  to  mind:  (i)  Is  there  an  optimal 
placement  of  fibers  in  a  matrix  and  (ii)  is  there  a  lower  limit  of  the  volumetric  fraction  of 
fibers  for  which  the  failure  of  a  composite  will  be  ideally  ductile  so  that  the  influence  of  the 
breakage  of  individual  fibers  may  be  ignored?  In  order  to  give  a  partial  answer  to  these 
important  questions,  we  investigated  the  following  problem. 

We  consider  a  cylindrical  fiber  which  is  embedded  into  a  matrix  binder.  Both  fiber 
and  matrix  materials  were  assumed  to  be  homogeneous  and  isotropic  but  of  different 
material  constants.  As  to  loading,  the  fiber  and  matrix  were  subjected  to  uniform  loads  ao 
and  oi  respectively  ( for  more  details  see  PART  VI).  The  analysis  revealed  the  load  transfer 
characteristics  from  the  matrix  to  the  fiber  and  vise  versa.  Theranalysis  was  subsequently 
extended  to  the  periodic  case  and  the  results  may  be  found  in  PART  VII.  Suppose  now  at  a 
certain  moment  in  time,  because  of  the  applied  loads,  a  plane  crack  perpendicular  to  the  axis 
of  the  cylinder  is  formed.  It  is  now  possible  to  estimate  th^  3D  stress  field  due  to  the 
presence  of  the  crack  (static  considerations)  and  thus  establish  the  zone  of  influence  of  a 
break  in  one  fiber.  These  results  can  then  be  extended  to  also  include  an  infinite  periodic , 
array  of  fibers.  The  results  are  presently  being  used  to  stuky  the  effect  of  briclging  inj 
cracked  composite  systems.  *  ~ 


5.  Adhesive  Butt  Joints. 


Finally,  the  problem  of  an  adhesive  butt  joint  has  been  investigated  and  the  results 
are  given  in  PART  Vin. 

6.  Gejsralizsd  Eailffl^-Critsria 


The  results  may  now  be  combined  in  order  to  derive  failure  criteria  for  crack 
initiation  applicable  to  unidirectional  composite  systems  based  on  3D  micromechanical 
considerations.  Much  of  this  work  has  already  been  carried  out  and  the  results  will  be 
reported  in  a  separate  report.  But  most  importantly  the  results  are  essential  for  a  better 
understanding  of  the  damage  evolution  process  in  critical  neighborhoods  such  as  holes, 
cracks,  joints,  edges  etc. 

For  example,  the  results  of  parts  I ,  II,  HI,  IV,  V,  VI,  and  VII  may  now  be  used  to 
address  the  phenomenon  of  bridging  in  ceramic  and  intermetalic  composites  and  thus 
develop  a  relationship  between  fracture  toughness,  fiber  volume  fraction  and  combined 
mechanical  thermal  loads. 
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PART  I 


ON  THE  THREE-DIMENSIONAL  STRESS  FIELD  OF  A  PERIODIC 

ARRAY  OF  FIBERS 
EMBEDDED  INTO  A  PLATE  MATRIX 


by 

E.S.  Folias  and  J.H.  Liu 


Department  of  Mathematics 
University  of  Utah 
Salt  Lake  City,  Utah  84112 


Abstract 

A  3D  micromechanical  model  has  been  developed  to  represent  a 
undirectional  composite  plate  which  is  subjected  to  a  uniform  trans¬ 
verse  load  (Tq.  The  model  assumes  the  fibers  to  be  cylindrical  inclu¬ 
sions  which  are  periodically  embedded  into  an  epoxy  matrix.  The 
materials  of  both  fibers  and  matrix  are  assumed  to  be  linear,  elastic, 
and  isotropic.  The  analytical  solution  shows  the  radial  stress  (x^r  to 
decrease  as  the  fiber  volume  fraction  Vj  increases.  The  stress  profile 
along  a  fiber  length  is  shown  to  be  constant  except  in  the  neighbor¬ 
hood  of  the  fiber  edge  where  a  boundary  layer  is  shown  to  prevail.  In 
this  re@on,  the  analytical  solution  shows  the  stress  field  to  be  singu¬ 
lar  which  is  a  departure  from  the  results  given  by  macromechanical 
theories. 

In  the  limit,  as  G/  —  0+  the  3D  stress  field  of  a  plate  weakened 
by  a  periodic  array  of  holes  is  recovered. 


1  INTRODUCTION 

It  is  well  recognized  that  fiber  composite  materials  are  very  attractive  for 
use  in  aerospace,  automotive  and  other  applications.  These  composites 
consist  of  relatively  stiff  fibers  which  are  embedded  into  a  lower  stiffness 
matrix.  Although  in  most  designs  the  fibers  are  aligned  so  that  they  are 
parallel  to  the  direction  of  the  external  loads,  it  is  almost  impossible  to 
avoid  induced  transverse  stresses  which  may  lead  to  premature  failure  of 
the  laminate.  An  excellent  example  of  this  is  the  case  of  a  filament  wound 
pressure  vessel  in  which  the  presence  of  curvature  induces  bending  as  well 
as  transverse  stresses  (Folias,  1965).  However,  in  order  to  be  able  to  predict 
their  failing  characteristics,  particularly  in  the  neighborhood  of  free  surfaces 
such  as  holes,  edges  etc.,  it  is  necessary  to  know  the  local  stress  behavior 
from  a  3D  point  of  view. 

An  overall  summary  of  some  of  the  results,  which  are  based  on  2D 
elasticity  considerations  can  be  found  in  the  books  by  Hull  (1981)  and  by 
Chamis  (1975).  In  their  pioneering  work,  Adams  and  Doner  (1967)  used 
finite  differences  to  solve  the  problem  of  a  doubly  periodic  array  of  elastic 
fibers  contained  in  an  elastic  matrix  and  subjected  to  a  transverse  load. 
Their  results  reveal  the  dependence  of  the  maximum  principal  stress  versus 
the  constituent  stiffness  ratio  {Ef  (Em)  for  various  fiber  volume  ratios.  A 
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few  years  later,  Yu  and  Sendeckyj  (1974)  used  a  complex  variable  approach 
to  solve  the  problem  of  multiple  inclusions  embedded  into  an  infinite  ma¬ 
trix.  Their  results  were  subsequently  specialized  to  cases  of  two  and  three 
inclusions  thus  providing  us  with  further  insight  into  the  strength  of  the 
composite. 

In  this  paper,  we  will  construct  an  analytical  solution  for  the  3D  stress 
field  of  a  matrix  which  has  been  reinforced  in  one  of  the  directions  with 
cylindrical  fibers. 

2  Formulation  of  the  problem 

Consider  the  equilibrium  of  a  body  which  occupies  the  space  |i|  <  oo,  \y\  < 
oo,  l^l  <  h  and  contains  a  periodic  array  of  cylindrical  inclusions  of  radius 
a,  whose  generators  axe  parallel  to  the  r-axis  (see  Fig.  1).  The  physical 
situation  depicted  here  is  that  of  a  uniderectional  composite  plate  that 
consists  of  a  matrix  where  fibers  are  embedded  into.  For  convenience,  all 
quanties  with  the  script  (m)  will  refer  to  the  matrix  while  quanties  the  with 
script  /  will  refer  to  the  fibers.  The  materials  of  both  matrix  and  fibers 
will  be  assumed  to  be  homogeneous,  isotropic  and  linear  elastic.  At  the 
interface,  i.e.  at  r  =  a,  perfect  bonding  will  be  assumed  to  prevail.  As 
to  loading,  a  uniform  tensile  stress  <To  is  applied  on  the  composite  plate 
(see  Fig.  1)  which  is  in  a  direction  perpendicular  to  the  axis  of  the  fibers. 
Furthermore,  the  surface  |r|  =  h,  for  both  regions,  matrix  and  inclusion, 
will  be  assumed  to  be  free  of  stress  and  constraints. 

In  the  absence  of  body  forces,  the  coupled  differential  equations  govern¬ 
ing  the  displacement  functions  are  : 

1  de^-'^  /,) 

+  ■  =  1-2.3;  >  =  ">,/  (1) 

where  is  the  Laplacian  operator,  Uj  is  Poisson’s  ratio,  u-"*^  and 
represent  the  displacement  functions  in  the  matrix  and  fibers  respectively, 
and 


e'«  = 


U) 


dxi 


i  =  1,2,3  ;  j  =  m,/. 


(2) 
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The  stress- displacement  relations  are  given  by  Hooke’s  law  as 


=  (3) 

where  Xj  and  Gj  are  the  Lame  constants  describing  the  material  properties 
of  the  matrix  and  of  the  inclusions. 

As  to  boimdary  conditions,  one  must  require  that  (see  Fig.  2  for  cell 
configuration) 

ai  \z\  =  h:  =  rjj)  =  0;  j  =  mj.  (4) 

at  r  =  a:  4"*)  -  =  ui*")  -  u</)  =  0.  (5) 

air)  -  -iP  =  -  r!P  =  -  r(p  =  0  (6) 

Moreover,  at  r  =  0  we  require  that  ail  stresses  and  displacements  be  fi¬ 
nite.  the  cell  configxuration  boundaries  AB  and  C D  will  be  taiken  as  planes 
of  symmetry,  thus  satisfying  the  respective  boundary  conditions  automati¬ 
cally.  It  remains,  therefore,  for  us  to  satisfy  only  the  continuity  boundary 


conditions  along  the  segment  BC,  i.e. 

„(■»>(«)  _  _  0)  + 

(7) 

-  f)  +  »'”•'(  j) 

(8) 

(9) 

(10) 

Fijmlly,  in  order  to  complete  the  formulation  of  the  problem, 
require  the  resultant  forces 

one  must  also 

/  Fxds  =  0 

Jbc 

(11) 
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/  Fyds  =  bag.  (12) 

Jbc 

It  is  found  convenient  at  this  stage  to  seek  the  solution  to  equation  (1)  in 
the  form; 

^  +  ;  i  =  1,2,3;  j  =  mj  (13) 

where  the  first  term  represents  the  particular  solution  and  the  second  term 
the  complementary  solution. 

3  Method  of  Solution 

A  generalized  ancilyticzil  solution  to  a  class  of  three-dimensional  problems 
which  arise  in  elastostatics  has  been  constructed  by  Folias  (1975  ,  1990a). 
The  solution  was  subsequently  used,  Folias  and  Wang  (1990b),  to  solve  for 
the  3D  stress  field  in  a  plate  which  has  been  weakened  by  a  hole.  This  work 
was  later  extended  to  the  solution  of  an  isotropic  inclusion  embedded  into 
a  matrix,  Penado  and  Folias  (1989).  On  the  bases  of  these  results,  one  may 
now  assume  the  solution  to  system  (1),  which  automatically  satisfies  the 
boundary  conditions  at  the  plate  faces  eq.  (4),  in  the  form  ^  : 


+A'')  - 

^  dx  mj  -h  1  dxdy 


^Note  that  because  of  symmetry  in  the  present  problem,  one  needs  only  to  consider  the 
region  0  <6  <  t|2. 
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(15) 


£  -d^ 


Zrrij  —  1  (j)  (j)  _  9^3  ^  _  1  ^ 

m^  + 1  ^  ^  5y  rrij  + 1  9x^ 


(16) 

From  which  the  stresses  can  easily  be  obtained  as; 


2^-12'''  =  ;;^  |{2/12^/.(A.) 

+  ^-3  -[2(mj  -  l)fi{0yz)  +  mjf2{0^z)]j 

00  r  dPmj)  ,  ,,  ,  , 

„?i  \~d?~ 

M!  _  ^^3^  2  1  2^^ 

^  5x  ^  5x*  rUj  +  1  5y  "*"  ruj  +  1  ^  dx^dy 

a3ffO)  ,dffb)  1 

-  (-^  -  ^^-^)l2(m,  -  l)/,(^.x)  +  m,/2(^.x)]| 

oo  ra^ffb)  jdHiiU  ,  , 

+  \~^  ■  '*’*~ar} 

2m,  aAji^^  _  ^  .  _  1  2^^ 

mj  +  1  5y  dx  ^  dx^  m^  +  1  ^  dx^dy 
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(19) 


J_^(C)0)  ^ _ ^  ~ 

2Gj  “ 


J  V 

mj-2  dx 


^Ml3.x) 


1  .(=)(,)  _  1  °° 


W/-‘  =  £  a?a?  +  m,/,(3.2)} 

oo  1  ,aHWl  ,  > 

'n?.  I  '  2“"~arj  ‘ 

,  m^  - 1  8a5'’  ^  a^  _  a‘4''  1 


where 


mj  + 1  dx  dx  dxdy  mj  + 1  dx^ 

2^^- =  -^2  I 

1  OO 


r(c)(i)  - 


TTlj  OO  a^ff.9) 


2Gj  mj  -  1  ax2 

1 00  a^^^J 
■^2  E 


{/3(;3„z)  +  /4(/?„2)} 
-  C03  (q„/i)  sin  (o„z), 


nsl 


HTT 


On  =  ~,n  =  1,2,3,. .. , 
h 


0^  are  the  roots  of  the  equation 


(20) 


(21) 


(22) 


(23) 


sin  {20^h)  =  -{20„h),  (24) 

.ffO)  ajid  ffU)  are  functions  of  x  and  y  which  satisfy  the  reduced  wave 
equation: 


a’  a*  ,  a^o) 

^ai2^ay2  ax 


(25) 
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3^  3^ 

dy 

(26) 

Aj^^  and  ^3^  are  two  dimensional  harmonic  functions,  and 

fii^^z)  =  cos  {$^h)  cos  {0^z) 

(27) 

f2{0z)  =  0uh  sin  {0uh)  cos  {0i>z)  —  0„z  cos  sin  (0^2) 

(28) 

f3{0uz)  =  cos  {0^h)  sin  {0^2) 

(29) 

f^{0u2)  =  0uh  sin  {0uh)  sin  {0^2)  +  0^2  cos  {0^,h)  cos  {0^z). 

(30) 

Examining  the  nature  of  the  boundary  conditions,  we  furthermore 
struct  the  solution  to  equations  (25)  and  (26)  in  the  form: 

con- 

=  g  akKk{0.r)  cos  (2k9) 

k=0 

(31) 

=  23  {okK2k(0yr)  +  ak/2*(^„r)}  cos  (2k0) 

kmO 

(31) 

»  r  cos  (2W) 

eSs 

(32) 

=  £  {cfciif3jk(a„r)  +  Ck/3fc(a„r)}  sin  {2k9) 

kmO 

(33) 

—  23  «»»  (2W). 

taO 

(34) 

(35) 

=  E  (4^  (2*  +  1)« 

kmO  ' '  * 

(36) 

E  +  C03  (2fc  +  i)« 

kmO  ' ^  ' 

(37)  i 

1 

1 

7 

(38) 


^  (-l)*Efcr  sin  {2k  +  1)9 

k=0 

=  g  (-l)*£fcr  C05  {2k  +  1)9  (39) 

k=0 

^  (— l)*Gfcr  cos  {2k  +  1)^  (40) 

At=0 

where  Ik{^r)  and  Ar*(/3r)represent  the  modified  Bess^  functions  of  the  first 
and  second  kind  and  ak,bic,Ck,dk.,  Ak,  Bk,  Ek  and  Gfc  are  arbitrary  constants 
to  be  determined  from  the  remaining  boirndwy  conmtions  (5)  -  (12).  Upon 
substitution  of  equations  (14)  -  (40)  into  equations  (5)  -  (12)  one  arrives 
at  a  system  of  twelve  equations  involving  series  in  z.  The  system  may 
then  be  solved  numerically  for  the  unknown  coefiicients.  Perhaps  it  may  be 
worth  noting  that  in  our  numerical  sinalysis  we  satisfied  first  the  boundary 
conditions  at  r  =  a  by  using  well  over  200  roots.  The  method  of  solution, 
as  well  as  the  rate  of  convergence  of  these  series,  is  similar  to  that  found  by 
Penado  and  Folias  (1989,  see  reference  for  details).  The  system  is  sensitive 
to  small  changes  and  for  this  reason  double  precision  was  used  throughout 
the  numerical  analysis. 
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4  Numerical  results 


Once  the  coefficients  have  numerically  been  determined,  the  stresses  and 
displacements  may  then  be  calculated  at  any  point  in  the  body.  Avoid¬ 
ing  the  long  and  tedious  numerical  details,  the  behavior  of  the  stresses  (t„ 
and  (Te$  versus  (z/h)  and  at  r  =  a  and  0  =  0  are  given  by  figures  3  and 
4  respectively.  It  is  noted  that  the  stresses  along  the  interior  length  of 
the  fiber  are  essentially  constant  and  that  as  one  approaches  the  edge  of 
the  fiber  length,  a  boundary  layer  is  shown  to  exist.  This  sudden  change 
suggests,  therefore,  the  presence  of  a  stress  singularity  at  such  regions.  In¬ 
deed,  a  separate  asymptotic  analysis  for  the  investigation  of  the  local  stress 
field,  at  such  neighborhoods,  shows  the  stresses  to  be  proportionzd  to  p~° , 
where  a  =  0.249  for  a  glass  fiber/epoxy  matrix  interface  (Folias  1989)  and 
Of  =  0.318  for  a  carbon  fiber  epoxy  matrix  interface  (Li  and  Folias  1990). 
It  may  also  be  noted  that  Figs  3  and  4  provide  us  with  important  infor¬ 
mation  concerning  the  regions  of  applicability  of  macromechanical  theories. 
The  reader  may  recall  that  such  theories  predict  the  stress  values  at  edges 
to  be  finite,  except  in  the  vicinity  of  an  interface  where  the  singularity 
strength  is  shown  to  be  very  weak  (Wang  and  Choi,  1982;  Folias,  1991). 
Thus,  they  tend  to  underestimate  the  actual  stress  levels  at  such  edges,  e.g. 
surface  of  a  hole,  surface  of  a  crack  etc.  But,  if  one  is  to  study  damage 
evolution  at  such  regions,  the  knowledge  of  the  local  stress  field  is  essen¬ 
tial.  Be  that  as  it  may,  a  closer  examination  of  Figs  3  and  4  shows  the 
boundary  layer  region,  for  a  transverse  applied  load,  to  be  restricted  to  a 
distance  of  one  fiber  diameter  away  from  the  fiber  edge.  On  the  other  hand, 
if  the  applied  lozul  is  in  the  direction  of  the  fiber  axis,  the  boundary  layer 
is  then  spread  out  to  a  distance  of  six  fiber  diameters  away  from  the  edge 
(Zhong  and  Folias  1991).  Thus  coupling  between  the  macromechanicsJ  and 
micromechanical  results  may  be  desirable  in  predicting  local  damage  due 
to  fracture. 

Returning  to  the  stress  profiles  (frr  and  (Tm  (Fig.  3  and  4)  we  note  that 
the  magnitude  of  the  stresses  decreases  as  the  fiber  volume  fraction,  V/, 
increases.  The  decrease,  however,  is  only  noticed  when  the  spacing  of  the 
fibers  becomes  less  than  four  fiber  diameters  center  to  center.  Figs  5  and 
6  show  typical  stress  profiles  for  Orr  and  a##  as  a  function  of  Gj/Gm-  It 

^In  this  analysis  the  maUrial  of  the  carbon  fiber  is  assumed  to  be  transversely  isotropic. 
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is  interesting  to  note  that  the  circumferencial  stress  osb  decreases  rapidly 
as  the  ratio  {Gj/Gm)  increases.  For  glass  fiber /epoxy  matrix  (G//Gm)  = 
16.67,  which  implies  that  the  (Tbs  stress  is  approximately  zero.  Thus,  all 
things  being  equal,  the  controlling  stress  for  possible  crack  failure  is  the 
radial  stress  tr,,  at  the  particular  location  of  ^  =  0.  For  large  {Gj/Gm) 
ratios,  the  radial  stress  reaches  an  asymptotic  value.  A  similar  result  was 
also  obtained  by  other  researchers  based  on  2D  considerations  (e.g.  Adams 
and  Donen  1967).  In  Fig.  7,  a  plot  d[  the  radial  stress  on  the  interface 
boundary  is  shown  to  decrease  as  the  angle  6  increases.  On  the  other 
hand,  the  shear  stress  TrB  (See  Fig.  8)  vanishes  at  the  two  positions  0  =  0“ 
and  90“  and  attains  its  maximum  value  at  0  %  45“ .  The  location  of  this 
maximum  shifts  slightly  to  the  right  as  the  ratio  of  (a/6)  increases.  Similar 
stress  profiles  also  appear  as  one  moves  towards  the  free  surface.  At  the 
free  surface,  the  question  arises  as  to  whether  the  strength  of  the  stress 
singularity  is  affected  as  the  separation  distance  between  fibers  becomes 
smaller  and  smaller.  While  initially  the  authors  believed  that  this  may  be 
the  case,  lately  they  believe  that  the  singularity  strength  will  not  be  altered 
but  that  the  function  multiplying  the  singular  term  is  expected  to  change. 
Be  that  as  it  may,  the  subject  is  under  further  investigation. 


5  Conclusions 

A  3D  micromechanical  model  has  been  developed  to  represent  the  response 
of  a  unidirectional  composite  plate  subject  to  a  transverse  load.  In  this 
model,  the  fibers  are  considered  to  be  cylindrical  inclusions  which  are  peri¬ 
odically  embedded  into  the  matrix.  The  material  of  both  fibers  and  matrix 
is  assumed  to  be  linear,  elastic  and  isotropic.  The  analysis  has  shown  that, 
as  the  fiber  volume  fraction  V/  increases,  the  radial  stress  (r^r  decreases  by 
30  to  40  percent.  On  the  other  hand,  the  circumferential  stress  om  is  almost 
negligible.  The  stress  profiles  across  the  fiber  length  are  almost  constant 
except  in  the  neighborhood  of  the  fiber  edge,  where  a  boundary  layer  is 
shown  to  prevail.  In  this  region,  the  stress  field  possesses  a  weak  stress 
singularity  which  for  a  glass  fiber/epoxy  matrix  composite  is  of  the  order 
0.25.  This  result  represents  a  departure  from  the  results  predicted  by  a 
macromechanical  theory.  This  inconsistency  is  attributed  to  the  fact  that 
macromechanical  theories  tend  to  average  the  local  effects  throughout  each 
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layer  thickness,  and  as  a  result  all  stresses  at  the  edge  are  predicted  to  be 
finite.  Thus,  the  present  analysis  also  provides  us  with  important  informa¬ 
tion  regarding  the  regions  of  applicability  of  macromechanical  theories. 

A  closer  examination  of  the  stress  field  reveals  that,  in  the  presence 
of  a  crack,  the  damaging  stress  for  possible  fzdlure  is  the  radial  stress 
particularly  at  the  location  ^  =  0.  Two  types  of  failure  immediately  come 
to  mind:  a  fiber/matrix  interface  crack  and  a  radial  matrix  crack.  In 
conjunction  with  this  work,  the  former  was  recently  considered  by  Folias 
(1991)  and  the  latter  by  Folias  and  Liu  (1991). 

It  is  well  recognized  that  void  nucleation  occurs  more  readily  in  a  triaxial 
tensile  stress  field,  a  result  which  is  consistent  with  experimental  observa¬ 
tions.  Such  a  model  for  estimating  the  void  nucleation  stress  may  now  He 
obtzuned,  if  in  our  previous  analysis  we  let  Gy  — >  O'*'.  The  physical  situation 
depicted  here  is  that  of  a  matrix  which  has  been  weakened  by  a  imiformly 
distributed  periodic  array  of  cylindrical  voids  or  holes.  Such  an  estimate 
may  then  serve  as  a  lower  bound  for  the  transverse  strength  of  a  unidirec¬ 
tional  composite  plate.  Without  going  into  the  nximerical  details,  we  plot 
in  Fig.  9  the  stress  concentration  factor  through  the  thickness  {z/h) 
and  for  a  typical  ratio  of  (a/6)  =  0.3.  It  is  noted  that  the  s.c.f.  is  relatively 
constant  throughout  the  interior  and  that  it  rises  slightly  as  it  reaches  the 
vicinity  of  the  free  surface  whereby  it  begins  to  drop  rather  abruptly.  The 
chzu’acteristic  stress  profile  is  in  agreement  with  that  found  by  Folias  and 
Wang  (1990)  for  the  case  of  one  hole,  the  variation  of  the  s.c.f.  as  a  fimc- 
tion  of  the  ratio  (a/6)  is  given  by  Fig.  10,  where  it  may  be  noted  that  the 
void  volume  fraction  in  the  matrix  is  given  by 

=  (41) 

It  is  clear  from  this  Fig.  that  the  s.c.f.  increases  rather  rapidly  as  the 
void  volume  fraction  ratio  increases.  The  result  is  in  agreement  with  our 
physical  expectations.  A  similar  stress  profile  is  also  observed  throughout 
the  thickness  including  the  plane  *  z  =  h,  (see  Fig.  11).  Finally,  in  Figs. 
12a  and  12b  we  plot  the  variation  of  the  s.c.f.  on  the  planes  z  =  0  and 
z  =  h  as  &  function  of  the  position  angle  $  and  for  different  (a/6)  ratios. 

^Folias  (1987)  has  shown  that  no  stress  singularity  is  present  in  the  vicinity  of  the 
intersection  of  the  hole  surface  and  the  free  of  stress  plane. 
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It  is  noted  that  in  the  region  70°  <  <  90“  the  s.c.f.  is  relatively  flat. 

This  is  the  region  where  a  crack  is  most  likely  to  initiate  and  subsequently 
propagate. 

In  closing,  it  may  be  appropriate  here  to  note  that  the  analysis  may  now 
be  extended  to  also  include  row  of  fibers  with  different  fiber  orientations, 
which  points  to  the  concept  of  a  laminated  composite  plate.  This  defines 
the  subject  of  a  subsequent  paper. 
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Fig.  1.  Geometrical  coofigmation 


Local  geometrical  configuration 
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Fig.  3.  .trra,  <7„  versus  thickness  (r/A)  »t  the  position 
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Fig.  5.  Radial  stress  versus  the  ratio  Gf/G^ 
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Fig.  6.  Circumferentizd  stress  versus  the  ratio  Gt/Gjn. 
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Fig.  7.  Radial  stress  versus  at  r  =  a. 
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Fig.  9.  The  stress  concentration  factor  versiis  the  thickness  {z/h)  for  a  peri¬ 
odic  2UTay  of  holes. 
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Fig.  10.  The  stress  concentration  {aM:tor 


versus  the  ratio  a/b,  at  z  =  0. 
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Fig.  12a.  The  stress  concentration  factor  versus  at  r  =  0. 
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Abstract 


This  paper  investigates  the  free  edge  effect  (xi  the  stress  field  of  a  carbon  fiber, 
which  is  embedded  into  an  epoxy  matrix.  The  fiber  is  assumed  to  possess  cylindrical 
symmetry  and  to  be  transversely  isotropic.  The  matrix  is  assumed  to  be  of  an  isotrr^ic 
material.  The  stress  field  is  induced  by  a  unifrxm  tension  applied  on  the  matrix  at  points  far 
away  from  the  fiber  surface. 

The  displacement  and  stress  fields  are  explicitly  derived  and  a  stress  singularity  is 
shown  to  prevail.  The  singularity  strength  is  shown  to  be  a  function  of  the  material 
constants  of  the  fiber  as  well  as  those  of  the  matrix.  Finally,  the  displacement  and  stress 
profrles  are  plotted  as  a  function  of  the  angle  <j>  which  is  measured  fitxn  the  free  surface. 


1.  Introduction 


Quite  often  in  engineering  practice,  structures  are  con^sed  of  two  elastic  materials  with 
different  properties  which  are  IxMided  together  over  some  surface.  Such  type  of  problem 
has  been  investigared  fiom  a  2D  point  of  view  by  many  researchers  and  the  results  can  be 
found  in  the  literature.  For  example,  Knein  (1927)  considered  the  plane  strain  problem  of 
an  orthogonal  elastic  wedge  btxided  to  a  rigid  base.  Rongved  (19SS)  investigated  the 
problem  of  two  bonded  elastic  half-spaces  subjected  to  a  ccxicentrated  force  in  the  interior. 
Subsequently,  ^^^ams  (1959)  studied  the  stress  field  arouixl  a  fault  or  a  crack  in 
dissimilar  media.  The  work  was  then  generalized  by  Rice  and  Sih  (1965)  to  also  include 
arbitrary  angles. 

It  was  not  until  1968  that  Bogy  (1968)  considered  the  general  problem  of  two 
bonded  quarter-planes  of  dissimilar  isotropic,  elastic  materials  subjected  to  arbitrary 
boundary  tractions.  The  problem  was  solved  by  an  application  of  the  Mellin  transfism  in 
conjunction  with  the  Airy  stress  function.  The  same  author  (1971)  extended  his  woric  to 
also  include  dissimilar  wedges  of  arbitrary  angles.  Shortly  thereafter,  Hein  and  Erdogan 
(1971),  using  the  same  method  of  solution,  independently  reproduced  the  results  by  Bogy. 
Finally,  Westmann  (1975)  studied  the  case  of  a  wedge  of  an  arbitrary  angle  which  was 
bonded  along  a  finite  length  to  a  half-space.  His  analysis  showed  the  presence  of  two 
singularities  close  to  each  other.  Thus,  elimination  of  the  first  singular  term  does  not  lead 
to  a  bounded  stress  field  since  the  second  singularity  is  still  present. 

Based  on  3D  consideration,  Luk  and  Keer  (1979)  investigated  the  stress  field  in  an 
elastic  half-space  containing  a  partially  embedded  axially-loaded,  rigid  cylindrical  rod.  The 
problem  was  formulated  in  terms  of  Hankel  integral  transforms  and  was  finally  cast  into  a 
system  of  coupled  singular  integral  equations  the  solution  of  which  was  sought 
numerically.  The  authors  were  able,  however,  to  extract  in  the  limit  from  the  integral 
equations  the  characteristic  equation  governing  the  singular  behavior  at  the  intersection  of 
the  free  surface  and  that  of  the  rigid  inclusion.  Their  result  was  in  agreement  with  that 
obtained  by  Williams  (1952)  fw  a  right-angle  comer  with  fixed-fiee  boundary  conditions. 

Haritos  and  Keer  (1979)  investigated  the  stress  field  in  a  half-space  containing  an 
embedded  rigid  block  under  conditions  of  plane  strain.  The  problem  was  formulated  by 
cleverly  superimposing  the  solutions  to  the  problem  of  horizontal  and  vertical  line 
inclusions  beneath  an  elastic  half-space.  By  isolating  the  pertinent  terms,  die  authors  were 
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able  to  extract  directly  from  the  integral  equations  the  oonder  of  the  stress  singularity  at  both 
comers.  Both  results  are  in  agreement  with  the  Willianos  solution.  Moreover,  the  authors 
point  out  the  importance  of  the  second  singularity  to  the  results  of  the  load  transfer 
problems. 

With  the  advancement  of  composites  and  their  extensive  use  in  the  aerospace 
industry,  problems  for  the  determination  of  the  stress  and  displacement  fields  around 
inclusions  have  drawn  considerable  attention.  For  instance,  Goodier  (1933)  investigated 
the  disturbing  effect  of  small  spherical  and  cylindrical  inclusions  (xi  an  otherwise  uniform 
stress  distribution  plate.  Numerical  results  were  presented  for  flaws,  perfectly-bonded 
rigid  inclusions  and  slag  globules  cases.  Hardiman  (1952)  used  the  complex  variable 
method  to  treat  the  elliptic  inclusion  problem,  he  was  the  first  to  find  that  a  unifOTm  applied 
load  at  infrnity  induces  a  ctxistant  state  of  stress  within  an  elliptic  inclusion.  The  work  was 
later  generalized  by  Sendeckyj  (1970)  to  include  the  solutions  of  the  elastic  curvilinear 
inclusitxi  problems.  It  was  not  until  1979  that  Tirosh,  Katz  and  Lifschuetz  (1979)  studied 
the  stress  interaction  of  a  single  frber,  embedded  in  an  elastic  matrix,  with  a  micro-crack 
situated  along  or  near  the  interface.  They  found  that  the  radial  tensile  stress  component  a„ 
is  higher  than  the  tangential  component  and  the  location  r  at  which  the  maximum  stress 
takes  place  is  not  on  the  interface  but  at  a  small  distance  ahead,  depending  (xr  the  Poisson's 
ratio  of  the  matrix.  Dundurs(1989)  noticed  that  the  stresses  in  a  body  that  contains  rigid 
inclusions  and  is  subjected  to  a  specified  surface  traction  depend  on  the  Poisson's  ratio  of 
the  materials.  If  the  Poisson's  ratio  is  set  equal  to  one  we  recover  the  case  of  plane  strain 
while  if  it  is  set  equal  to  infinity  we  recover  the  case  of  plane  stress. 

In  the  area  of  multiple  or  periodical  inclusions  problems,  Adams  and  Doner  (1967) 
obtained  a  2D  nunKiical  solution  by  a  systematic  overrelaxation  procedure  for  a  plate 
containing  a  rectangular  array  of  inclusions  embedded  in  an  elastic  matrix  and  subjected  to  a 
uniform  normal  transverse  stress  at  infrnity.  At  the  same  time,  Goree  (1967)  presented  a 
solution  for  the  stress  and  displacement  distributions  in  an  infinite  elastic  matrix  ccxitaining 
two  perfectly-bonded  rigid  cylindrical  inclusions  of  different  radii.  Subsequently,  Haeno- 
and  Ashbaugh  (1967)  used  the  displacement  potentials  method  to  express  die  stress 
distributicxis  in  a  unidirectional  multiple  fibers  composite  plate  under  external  and  residual 
loads.  Marloff  and  Daniel  (1969)  used  a  standard  stress-freezing  technique  to  determine 
the  3D  stress  distribution  in  the  matrix  of  a  unidirectional  composite  plate  subjected  to 
matrix  shrinkage  and  normal  transverse  load.  Based  txi  Sendeckj's  previous  work,  Yu  and 
Sendeckj  (1974)  extended  their  study  by  means  of  the  Schwarz  alternating  method  to  the 
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case  of  an  iniinite  elastic  matrix  containing  random  number  of  elastic  inclusitms. 

Moreover,  Adams  and  Crane  (1984)  studied  a  microscopic  region  of  a  unidirectional 
composite  plate  by  the  finite  element  micromechanical  analysis  using  a  generalized  plane 
strain  formulation.  Finally,  Keer,  Dundurs  and  Kiattikomol  (1973)  studied  the  phenomena 
of  the  separation  of  a  smooth  circular  inclusicm  fiom  a  matrix  which  is  subjected  to  an 
uniform  load.  Using  finite  integral  transforms,  die  problem  of  finding  the  extensicm  of 
separation  and  the  contact  pressure  is  reduced  to  the  solutitxi  of  a  Fredholm  integral 
equation  with  a  weakly  singular  Kernel. 


Recently,  Folias  (1989)  examined  the  local  stress  field  in  the  neighborhood  where  a 
fiber  embeddded  into  an  epoxy  matrix  meets  a  fiee  surface.  In  this  analysis,  Folias 
assumed  the  fiber  and  the  matrix  to  be  isotropic  but  of  different  material  constants.  The 
analysis  showed  that  the  stress  field  in  this  vicinity  to  be  singular.  Moreover,  the  3D 
analysis  showed  that  to  the  case  of  a  cylindrical  inclusion,  tl^  order  of  the  singularity 
strength  is  precisely  that  which  was  reported  by  Bogy  (1968)  based  (xi  2D  consideraticxi. 
In  the  present  paper  the  analysis  is  extended  to  also  include  a  transversely  isotropic  fiber, 
e.g.  carbon  fiber,  and  the  explicit  3D  displacement  and  stress  fields  are  recovered. 

2.  Formulation  of  the  problem 


Consider  the  equilibrium  of  a  cylindrical  carbon  fiber  of  radius  a  which  is 
embedded  into  a  hoirx>geneous,  isotrc^ic  and  linea’  elastic  matrix  that  occupies  the  space 
|x|<oo,  |yl<oo  and  Ul^h.  The  fiber  is  assumed  to  be  of  a  transversely  isotropic 
material  with  different  material  properties  than  those  of  the  matrix.  Moreover,  tire  axis  of 
the  fiber  is  assumed  to  intersect  the  bounding  plane  z^  ±h  perpendicularly  and  that  the 
matrix  is  subjected  to  a  urtiform  tensile  load  Go  in  the  direction  of  the  y>axis  and  parallel  to 
the  xy-plane  (see  Fig.  1).  Perfect  bonding  at  the  interface  is  assumed  to  prevail. 


In  the  absence  of  body  forces,  the  equilibrium  equations,  in  terms  of  the  stresses 
Gjj,  are: 


dt  X  de  di  X 

Bi  X  do  Bz  X 


Bx  X  BO  Bz 


=  0 


(1) 

(2) 

(3) 
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where  superscript  m=l  represents  the  material  of  the  matrix  and  m=2  represents  the 
material  of  the  fiber. 

The  stress-strain  relations  for  a  transversely  isotropic* ,  as  well  as  an  isotropic, 
material  are  given  by  the  ctmstitutive  relations 
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where  the  represent  the  respective  material  constants. 

As  to  boundary  conditions,  we  require  that  (i)  the  stresses  on  the  planes  ±h  must 
vanish  and  (ii)  that  the  displacements  and  stresses  of  the  two  materials  must  match  at  the 
interface,  i.e.  at  p=a. 


3.  Method  of  solution 


The  primary  objective  of  this  analysis  is  to  derive  an  asymptotic  solution  which  is 
valid  in  the  immediate  vicinity  of  the  comer  points,  i.e.  the  points  where  the  interface  meets 
the  free  surface  of  the  matrix.  Thus  following  the  same  method  of  solutitm  as  that  of  Folias 
(1989),  one  may  write  the  equilibrium  equations  (l)-(3)  in  terms  of  the  displacements  Ujjf™) 
as^ 


r  <■>  I  (r  <-i  +  c  ic  1-1  ^“,"’-0 


( 


Si — l15s_+ —  0 


Ja(r-a)* 


dz^ 


C J”*  +  (€44^"*  +  + C, 


^(r-a)^ 


^(r-a)^ 


'33 


=0 


(5) 

(6) 

(7) 


*  The  material  of  the  fiber  is  assumed  to  possess  a  cylindrical  symmetry. 
2  In  writing  the  following  equations,  we  assumed  that  (r-a) «  a. 
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It  is  interesting  to  note  that  die  stress  field,  at  the  vicinity  of  the  fiber  interface,  leads 
to  two  coupled  equations  (5)  and  (7)  for  the  displacements  and  w<™)  and  an  additicnial 
equation  (6)  for  the  displacement  ue(*n). 


Without  going  into  the  mathematical  details,  by  direct  substitution,  me  can  show 
that  the  following  displacement  field  satisfies  die  governing  equations  (S)-(7): 


d(T-a)dz 


(«n)  ^ _ ./-I  (m)  ^ 

"  "r"  *=r 


u^(«)  ^  H(«) 


where  the  psudo-harmonic  functions  H*"*  and  satisfy  the  equations: 
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(13) 


(14) 


(15) 


(16) 


It  is  found  convenient  at  this  stage  to  introduce  the  local  coordinate  system 
(see  Fig.  2), 

r-a  =  pcos^ 
h-zspsin^ 

and  to  adopt  the  following  definitions: 


(17) 

(18) 
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(28) 

(29) 

8 


Substiniting  the  previously  ccmstructed  displacement  field  into  die  boundary 
conditions  (26)-(29),  we  arrive  at  a  system  of  twelve  algebraic  equations,  the  determinant 
of  which  must  vanish.  This  latter  condition  leads  m  a  transcendental  equation  for  die 
characterisdc  value  a. 

Without  going  into  the  mathematical  details,  these  characteristic  values  a  may  easily 
be  determined  with  the  aid  of  a  computer.  Althou^  the  transcendental  equation  has  an 
infinite  number  of  real  toots,  only  those  values  which  lie  in  the  interval  1<  a  <  2  are  of 
practical  interest  Furthermore,  in  this  interval  diere  are  no  complex  roots  present  In 
general  a  depends  on  the  respective  material  prc^ierties  of  the  matrix  as  well  as  of  the  fiber. 

The  displacement  and  stress  fields  may  now  be  computed  in  terms  of  one  unknown 
constant  which  in  turn  is  to  be  determined  fiom  the  loading  conditions  far  away  fixxn  the 
fiber-matrix  interface.  Without  going  into  the  mathematical  details,  the  explicit 
displacement  and  stress  fields  are  found  to  be: 


(i)  displacement  field: 


'  -r4r(c..‘"’ + Cu'-')a(a + cosifiB) 


■  -A,<->siii[(o-lW,]+B,<->cos[(a-l)*,]+--^|  v^">(0 

V  /  0 

cos[(a  -  l)0j-(a  +  l)C]dC  - 

a(a+l)  J  (30) 


w‘“’  =  a{a  +  1)P2‘'"‘  cos(^0) 


■  [c..'"’  -  ^]{ A,‘-‘cc«{(o  -  !)♦,] + B  “  »m[(o  -!)♦,] 


(31) 
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=  o{a  +  cos(fi0) 

cos[(a  - 1)^,] + sin[(a  - 1)^.]} + 0(a“) 


(32) 


(ii)  stress  Held: 

=  — J— C44‘”’(a  -  l)a(of  +  cosifiB) 

sta[(a  -  2W.]  -  B,‘-'  cos[(a  -  2)#,] 
i— I  *<<•>  ({:)cos[(o  -  2)^  -  (a  +  Dfjiin  +  [(C.," 


(a  +  1) 


C,/"^  sin20,cos^a  ^  ,,,  sin>,  ^  cos^^.sin^J 

^  a(a  +  l)  ”  “  fle(a  +  l)  £i‘”^  a(a  +  l)  J 
^  [  "  (a-l)a(a  +  l)  (a-l)a(a  +  l)J 

=  — rW(“  "  DPj*'*  COSO50) 


C,  «(C, +  c«<-’)  -  c.,'->fc.,'->  -  ^ j  {A,'"> 
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cos[(a  -  2W,  -  (a + DfK} + {C..“’(C„‘-' + C."”)  -  C„ 

(mi  ^44  sin2^2®®*^2  +  r  r  ("Jfp  ^“'  +  0 

I  ”  "^)\  a(a  +  l)  ^  " 


Ja(a+1)  a(a+l)  T 
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(33) 


(34) 
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co&(fi9) 

I  C„<->(C„‘->  +  C«'*>)-C„<->[c.,->-^j  {A,'-’ 

sin[(o  -  2)*,]  -  B  “  cos[(o  -  2W.1-  I  V^”(f  > 

cos[(o  -  -  (o  +  1)CK}  +  {[C„“’(C.,'-’  +  C  J->)  -  C„<"> 

(c  <“)  -  £is!ril  sin  2^2  cos  0;  ^  [-c,  <“>  ((:,,*“'  +  )  + 

I  “  JJ  a(a  +  l)  ^ 

(-)c  (,)i 

"  ^a(a  +  l)  a(a  +  l) 
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f  <“>  (a-DaCa  +  Dj"^ 


=  €44^“’ (a  -  l)a(a + Dpj""^  cos^d) 


.  j{A2‘“’  cos[(a  -  2)0j] + Bj<"’  sin[(a  -  2)0^] 

+ — - —  I  \^^"’(Osin[(a  -  2)02  -  («  +  l)C]dc| 
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t =  — pi— rC44^“’(a  - 1)0(0 + l)p"'^cos(fi0) 

*  -  /s  (“) 


{-A,^”’  sin[(o  - 1)^,]  +  B/”’  cos[(o  - 1)^,]} + 0(p,®‘‘ ) 


(35) 


(36) 


(37) 


(38) 


=  9iLll£xi!l(^a  -  l)a(a  +  l)p  **cosO30) 
{A,<“>cos[(a  -  l)0j + B»<">sin[(a  - 1)0,]}  +  OCPi*’ ) 


4.  Discussion 

As  a  practical  matter,  two  different  types  of  fibers  embedded  into  an  epoxy  matrix 
wens  considered,  a  carbon  fiber  and  a  glass  fiber  defined  by  the  following  elastic 
properties  (see  Table  1.): 

Table  1.  The  Material  Properties  for  Carbon  Fiber,  Glass  Fiber  and  Epoxy  Matrix 

Carbon  Fiber  Glass  Fiber  EP0?ty  Mam 

Cii=20.40MPa  Cii=  99.19  MPa  Cii=  6.62  MPa 

Ci2=  9.40  MPa  Ci2=  27.69  MPa  Ci2=  3.4 1  MPa 

Ci3=10.50  MPa 
C33«240.00  MPa 
C44S  24.(X)  MPa 

Omitting  the  long  and  tedious  numeiical  details,  the  characteristic  values  of  a  were 
found  to  be  0=1.693  and  0=1.737,  respectively^.  The  analysis  suggests  that  the  presence 
of  a  carbon  fiber  induces  a  slightly  higher  singular  stress  field  than  that  of  a  glass  fiber  and 
consequently  is  more  prone  to  failure.  Moreover,  in  the  limit,  one  recovers  precisely  the 
corresponding  isotropic  solution  derived  by  Folias  (1989).  In  Figs.  3-6,  the  behaviors  of 
the  local  displacement  and  stress  fields  as  a  function  of  angular  distribution  0  are  depicted 
for  both  fibers  considered  above.  The  reader  may  note  that  for  simplicity  we  have  adopted 
the  definition: 

C{9)  =  A2”’a(a+l)cos(/J0)  (39) 

which  in  the  limit,  as  the  material  constants 

Cg  .  . 

isotropic,  (40) 


3  It  may  be  noted  that  a  second  real  root  exists  widiin  this  interval,  howevo*  it  leads  to  a 
slightly  weaker  stress  singularity  (see  Westmann,  1975). 
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Rg.  4  Displacements  versus  ^  for  transversely  isotropic  case 
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stresses 


stresses 


the  function 


0(6)  -4  C*(0) 


(41) 


It  is  noteworthy  to  note  that,  although  the  displacement  and  stress  profiles  correspcmding  to 
the  two  different  fibers  are  very  similar,  their  respective  magnitudes  are  different  Similar 
trends  are  also  expected  to  prevail  at  locations  further  away  from  the  free  surface.  Finally, 
it  remains  for  us  to  relate  the  unknown  ctmstant  to  the  applied  load  far  away  from  the 
vicinity  of  the  fiber.  This  matter  is  presently  under  investigation,  and  the  results  will  be 
reported  in  a  future  paper  where  we  will  also  take  into  account  other  effects  such  as  stresses 
due  to  a  temperature  mismatch,  mechanical  loads  along  and  perpendicular  to  the  fiber  axis, 
and  a  correction  factor  to  account  for  the  presence  of  a  periodic  array  of  fibers  embedded 
into  a  matrix.  Finally,  the  variation  of  the  exponent  a  as  a  function  of  the  material  ratio 
033(2)7033(1)  is  ^ven  by  Fig.  7. 
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ABSTRACT 


This  paper  deals  with  the  3D  stress  field  of  a  cylindrical  fiber  which  is 
embedded  into  a  resin  matrix.  The  composite  is  then  subjected  to  a  uniform 
tensile  load  oq.  The  strain  energy  retease  rate  is  computed  and  the  criterion  is 
used  to  predict  debonding  initiation  at  the  fiber/matrix  interface.  The 
analysis  shows  that  this  fidlure  is  most  likely  to  occur  at  the  free  surface,  ie  the 
region  where  the  fiber  intersects  a  free  surface  for  example  a  hole,  an  edge,  or 
a  crack.  Moreover,  it  will  occur  at  approximately  (1/10)  the  load  value 
required  for  the  same  failure  to  commense  at  the  center  of  the  fiber  length. 

The  results  are  also  extended  to  include  a  doubly  periodic  array  of  fibers 
which  are  embedded  into  a  matrix.  Based  on  3D  considerations,  the  stiffness 
matrix  is  shown  to  increase  as  the  volume  fraction  of  the  fibers  increases. 
Similarly,  the  stress  Orr  in  the  matrix  is  shown  to  decrease  as  the  volume 
fraction  of  the  fibers  inaeases. 
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INTRODUCTION 


It  is  well  recognized  that  fiber  composite  materials  are  very  attractive 
for  use  in  aerospace,  automotive  and  other  applications.  These  composites 
consist  of  relatively  stiff  fibers  which  are  embedded  into  a  lower  stiffness 
matrix.  Although  in  most  designs  the  fibers  are  aligned  so  that  they  are 
parallel  to  the  direction  of  the  external  loads,  it  is  almost  impossible  to  avoid 
induced  transverse  stresses  which  may  lead  to  premature  failure  of  the 
laminate.  An  excellent  example  of  this  is  the  case  of  a  filament  wound 
pressure  vessel  in  which  the  presence  of  a  curvature  induces  bending  as  well 
as  transverse  stresses  (Folias,  1965).  However,  in  order  to  be  able  to  predict 
their  failing  characteristics,  particularly  in  the  neighborhood  of  free  surfaces 
such  as  holes,  edges  etc.,  it  is  necessary  to  know  the  local  stress  behavior  from 
a  3D  point  of  view. 

An  overall  summary  of  some  of  the  results,  which  are  based  on  2D 
elasticity  considerations  can  be  found  in  the  books  by  Hull  (1981)  and  by 
Chamis  (1975).  In  their  pioneering  work,  Adams  and  Doner  (1967)  used  finite 
differences  to  solve  the  problem  of  a  doubly  periodic  array  of  elastic  fibers 
contained  in  an  elastic  matrix  and  subjected  to  a  transverse  load.  Their 
results  reveal  the  dependence  of  the  maximum  principal  stress  versus  the 
constituent  stiffness  ratio  (Ef/Em)  for  various  fiber  volume  ratios.  A  few 
years  later,  Yu  and  Sendeckyj  (1974)  used  a  complex  variable  approach  to 
solve  the  problem  of  multiple  inclusions  embedded  into  an  infinite  matrix. 
Their  results  were  subsequently  specialized  to  cases  of  two  and  three 
inclusions  thus  providing  us  with  futher  insight  into  the  strength  of  the 
composite.  On  the  other  hand,  the  separation  of  a  smooth  circular  inclusion 
from  a  matrix  was  investigated  by  Keer,  Dimdurs  and  Kiattikomol  (1973).  By 
using  finite  integral  transforms,  they  were  able  to  reduce  the  problem  to  that 
of  a  Fredholm  integral  equation  with  a  weakly  singular  kernel.  Thus, 
extracting  the  singular  part  of  the  solution,  they  were  able  to  reduce  the 
remaining  problem  to  a  simpler  one  which  lends  itself  to  an  effective 
numerical  solution.  Their  results  are  very  general  and  are  applicable  to 
various  combinations  of  material  properties  and  loads. 
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In  this  paper,  use  of  the  local,  3D,  stress  field  will  be  made  in  order  to 
examine  the  dependence  of  the  stress  Crr,  in  the  matrix,  on  the  ratio  (Gf/Gm). 
The  strain  energy  release  rate  will  then  be  computed  in  order  to  predict  crack 
initiation  at  the  fiber/matrix  interface.  Particular  emphasis  will  be  placed  in 
the  region  where  fibers  meet  a  free  surface  as  well  as  at  the  center  of  a  fiber's 
length. 


FORMULATION  OF  THE  PROBLEM 


Let  us  consider  a  cylindrical  fiber  of  homogeneous  and  isotropic 
material,  e.g.  a  glass  fiber,  which  is  embedded  into  a  matrix  of  also 
homogeneous  and  isotropic  material. 

Futhermore,  we  assume  the  matrix  to  be  a  rectangular  plate  with  finite 
dimensions  2w,  2/,  and  2h  as  defined  by  fig.  1.  For  simplicity,  we  assume 
w  / 

—  >  8  and  ^  >  8.  Such  an  assumption  will  guarantee  that  the  boundary  planes 

X  =  ±  w,  and  y  =  ±  /,  will  not  effect  the  local  stress  field  adjacent  to  the  fiber.* 
Thus,  mathematically,  one  may  consider  the  boundaries  in  the  x  and  y 
directions  to  extend  to  infinity.  As  to  loading,  the  plate  is  subjected  to  a 
uniform  tensile  load  Oq  in  the  direction  of  the  y-axis  and  parallel  to  the 
bounding  planes  (see  Fig.  1). 

In  the  absence  of  body  forces,  the  coupled  differential  equations 
governing  the  displacement  functions  u^i^  are 


1 

1-2  vj 


ae<j) 

dxi 


+  V2 


=  0,  i  =  1,2,3,  j  =  1,2, 


(1) 


where  V2  is  the  Lapladan  operator,  Vj  is  Poisson's  ratio,  up^  and  up^  represent 
the  displacement  functions  in  media  1  (matrix)  and  2  (fiber)  respectively,  and 


e(j) 


1,2,3;  j  =  l,2. 


(2) 


The  stress-displacement  relations  are  given  by  Hooke's  law  as 


*  This  can  be  seen  from  the  results  which  were  recently  reported  by  (Penado  and  Folias  (1989). 
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where  Xj  and  Gj  are  the  Lame'  constants  describing  media  1  and  2. 

THE  SOLUTION  FOR  ONE  HBER 

A.  Region  where  fiber  intersects  the  free  edge 

This  problem  was  recently  investigated  by  the  author  (Folias  1989)  who 
was  able  to  recover,  explicitly,  the  three  dimensional  stress  field  adjacent  to 
the  surface  of  the  fiber*  .  Without  going  into  the  mathematical  details,  the 
displacement  and  stress  fields  for  the  matrix  are  given  in  terms  of  the  local 
coordinate  system  (see  fig.  2)  by: 

(i)  displacement  field: 

u(^)  =  An  sine  {b  [2  (1-vi)  cos  (a-l)(|>  -  (a-1)  sin^  sin(a-2)())]  (4) 

-  (a+1)  [(l-2vi)  sin  (a-1)^  +  (a-1)  sin0  cos  (a-2)  ()>] } cos(2n0) 

v(l)  a  An  COS0  {b  [2  (1-vi)  cos  (a-l)(l)  -  (a-1)  sin(|>  sin(a-2)(t)  ]  (5) 

*  (a+1)  [(l-2vi)  sin  (a-1)  ^  +  (a-1)  sin(|>  cos  (a-2)0]  }  cos(2ne) 

w^^^  a  An  p“"^  {b  [-  (l-2vi)  sin  (a-l)<>  +  (a-1)  sin4)  cos  (a-2)  <J)  ]  (6) 

-  (a+1)  [  2  (1-vi)  cos  (a-l)4>  +  (a-1)  sin0  sin(a+2)  0  ]  }  cos(2ne) 

(ii)  stress  field: 

2G^^)  (a-1)  An  p®"^  (b  [  2  cos  (a-2)(J>  -  (a-2)  sin(t»  sin(a-3)(t)]  (7) 

-  (a+1)  [sin  (a-2)(|>  +  (a-2)  sin0  cos  (a-3)^  ]  }  cos(2ne) 

-  4viG(^^  (a-1)  An  p®’^  {b  cos  (a-2)0  -  (a+1)  sin(a-2)(|>}  cos(2n0)  (8) 

*  A  similar  analysis  for  a  transversely  isotropic  fiber  meeting  a  free  surface  has  recently  been 
completed  and  the  results  will  be  reported  soon. 
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(9) 


2G(1)  (a-1)  An  p«-2  {b  (a-2)  sm^  sm(a-3)(|> 

+  (a+1)  [(a-2)  sin^  cos  (a-3)^  -  sin(a-2)(>  ]}  cos(2n0) 

*  2G^^)  (a-1)  An  p“‘2  (b  [sin  ia-2)^  +  (a-2)  sin(|>  cos  (a-3)(t>] 


(10) 


-  (a+1)  (a-2)  sin<>  sin(a-3)4> 


}cos(2n0) 


(1)  (1)  „ 
^r0  =  '^02  =  0' 


(11)-(12) 


where  n  =  0,1/2,...  and  B  is  a  function  of  the  material  constants  and  An  is  a 
constant  to  be  determined  from  the  boundary  conditions  far  away  from  the 
fiber*.  In  general,  the  characteristic  value  of  a  depends  on  the  material 
constants  of  the  fiber  as  well  as  of  the  matrix.  A  typical  example  is  given  in 
fig.  3. 

Upon  examination  of  the  stress  field,  the  following  remarks  are  worthy 
of  note.  First,  the  stress  field  in  the  neighborhood  where  the  fiber  meets  the 
free  surface  is  signular.  Moreover,  in  the  limiting  case  of  a  perfectly  rigid 
inclusion  this  singularity  strength  reaches  the  value  of  0.2888.  Second, 
boundary  conditions  and  satisfied  as  a  consequence  of  the 

odd  functional  behavior  in  <|>,  which  points  to  the  presence  of  a  boundary 
layer  solution  as  one  approaches  the  free  surface.  Third,  on  the  free  surface 
the  radial  stress  is  (1/vi)  times  the  drcumpherential  stress.  This  suggests, 
therefore,  that  if  a  crack  was  to  initiate,  it  would  propagate  along,  (or  very 
very  close  to)  the  fiber/matrix  interface.  Clearly,  the  occurance  of  either 
adhesive  or  cohesive  failure  will  depend  on  the  relative  strengths  of  the 
interface,  of  the  fiber,  and  of  the  matrix.  All  things  being  equal,  the  analysis 
shows  the  stresses  to  be  highest  at  the  interface,  thus  pointing  to  an  adhesive 
type  of  failure. 


*  for  one  fiber  n»0,l.,  while  for  a  periodic  extension  n»0,l,2,.. 
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6.  Interior  region 

The,  3D,  stress  field  for  this  region  has  also  been  recovered  by  Penado 
and  Folias  (1989)  and  the  results  for  various  (a/h)  and  (G2/G1)  ratios  may  be 
found  in  the  literature.  The  results  have  subsequently  been  extended  (Folias 
and  Liu,  1990)  to  also  include  a  layer  of  modified  matrix  around  the  fiber. 
Thus  for  VI  =  0.34,  V2  =  0.22  and  (G2/G1)  =  16.67  the  stresses  a^^^and  at  r  = 

a  and  for  all  I  z  I  2  h  are  given  in  figs  4  and  5  respectively.  Finally,  fig.  6  (for 
X.=0)  shows  the  variation  of  the  stress  o^^as  a  function  of  the  ratio  (G2/G1). 

INTERFACE  FAILURE  CLOSE  TO  THE  FREE  SURFACE 

A  closer  inspection  of  the  local  stress  field  shows  that  a  crack  is  most 
likely  to  initiate  at  the  location  6  =  0  and  subsequently  propagate  along  the 
fiber/matrix  interface  until  it  reaches  a  nominal  value  of  the  arc  length 
beyond  which  it  will  advance  into  the  matrix.  Moreover,  once  the  crack 
begins  to  propagate,  it  will  simultaneously  propagate  along  the  interface  and 
parallel  to  the  axis  of  the  fiber  (mode  HI).  Thus,  crack  propagation  will  be 
governed  initially  by  a  mode  I  failure  and  subsequently  by  a  combination  of 
mode  I  and  mode  III  failure.  It  is  now  possible  for  us  to  examine  the  first 
stage  of  the  failing  process  and  to  obtain  an  estimate  of  the  debonded  arc 
length  as  well  as  an  estimate  of  the  critical  transverse  stress  for  crack 
initiation. 

As  a  practical  matter,  we  will  consider  the  special  case  of  a  glass  fiber 
embedded  into  an  epoxy  matrix  with  the  following  properties 

Gi=2.10GPa  VI  =0.34  (13) 

G2  =  35.00  GPa  V2  =  0.22. 

Without  going  into  the  numerical  details,  the  constants  a,  A  and  B  for  this 
example  are  found  to  be’ : 

*  The  constant  A  has  been  determined  by  comparing  the  displacement  ^  as  well  as  the 

stress  at  0  =  0,  at  z  s  h  and  for  (a/h)  =  0.5  with  the  work  of  Penado  and  Folias  (1989). 

A*  Ia„. 
n 
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a  =  1.7511,  /i  aa-2  -  0.6349  Oo ,  B  -  2.1302  , 


(14) 


where  Oq  now  has  the  units  of  GPa.  Thus,  from  equations  (7)*(12)  one  has 
(i)  at  <>  =  0  and  6  =  0: 


-  4.0633  oo 


0.2489 


.(1) 


.(1) 


(ii)  at  (|>  =  n/2  and  0  =  0: 

Jl) 


rr 

’ee 

rz 


1.9163  00  f^JO.2489 


0.4844  o. 


(1) 


-  0.2930  o 


rr 

(1) 


rr 


(15) 

(16) 

(17) 

(18) 
(19) 


It  is  dear  now  from  equations  (15)  and  (16)  that  crack  failure  is  most 
likely  to  initiate  and  subsequently  propagate  along  the  fiber/matrix  interface 
rather  than  perpendicular  to  it.  Similarly,  equations  (17)  and  (19)  suggest  that 
failure  in  the  direction  parallel  to  the  axis  of  the  fiber  is  dominated  first  by  a 
mode  I  and  second  by  a  mode  III  type  of  failure.  It  may  also  be  noted  that 
Oyy^  attains  a  maximum  at  6  =  0  and  descreases  as  one  travels  along  the 
surface  of  the  fiber. 

Finally,  based  on  3D  considerations,  the  stress  field  away  from  the 
edges,  2  =  ±h,  and  in  the  interior  of  the  plate  was  shown  to  be  non-singular 
(Penado  and  Folias  1989,  Folias  and  Liu  1990)  with* 

=  0.4090  =  0.4090  (1.4281  Co)  =  0.5841  Oo ,  at  6  =  0  (20a) 


•  W 

These  results  are  valid  for  a  ratio  of  (a/h)  «  0.05  and  subject  to  the  assumption  that  ("J )  >  8 

and  (//a)  >  8  in  which  case  the  end  boundaries  in  x  and  y  have  insignificant  effects  on  the 
local  to  the  fiber  stress  field. 
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at  r  =  a  and  2  =  0.  Comparing  this  value  with  that  of  the  corresponding  plane 
strain  solution 


(1)  VI  (1) 


(20b) 


one  notices  that  it  is  approximately  13%  higher  in  value  due  to  the  presence 
of  the  stresses  in  the  third  dimension. 

It  is  now  possible  for  us  to  obtain  an  approximate  criterion  for  de¬ 
bonding  along  the  fiber/matrix  interface  based  on  Griffith's  theory  of  fracture. 
Thus,  following  the  work  of  Toya  (1974),  if  one  assumes  the  presence  of  an 
interface  crack  of  length  2aP  and  if  futhermore  takes  into  account  the  local  3D 
stress  field,  then  Toya's  result  may  be  written  as 

(1/16)  (1.1337ao  )2  k  a  Ai  (1+4  e2)  n  No  Nosin  (5  exp  [2  e  (ti-P)]  =  2y  (21) 


where 


l+k2 

l+k2  +  (l+k])  (G2/G]) 


3-4vi  for  plane  strain 

ki  =  3  -  Vi 

^  for  plane  stress 


Ai  = 


{i+k]  + 


a+k2)  (Gi/G2) 


=  "27c  Lki+C 


1  +k2  (G1/G2) 
kl  +(Gi/G2) 


KT  1  ^  l+k2(Gi/G2)  r  , 

No  =  Go‘j^-  k  ki  +  (Gi/G2) 


Go  = 


1-  (cos  p  +  2£  sin  p)  exp  [2£  (it-p)]  +  (1-k)  (1  +  4e2)  sin^p 
2-k-k  (cos  P  +  2£  sinp)  exp  [2£  (n  -  P)] 
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where  No  is  the  complex  conjugate  of  No,  "iii  is  the  specific  surface  energy  of 

the  interface  and  P  the  angle  of  debonded  interface  (see  Fig.  7).  While  it  is 
true  that  this  type  of  approach  does  not  provide  results  for  the  exact  initiation 
of  an  interface  crack  problem,  ie  from  a  condition  of  perfectly  bonded  interface 
to  that  of  a  partially  debonded  interface,  it  does,  however,  provide  a  very  good 
first  approximation  to  this  complex  phenomenon.  The  author  is  well  aware 
of  that  and  is  presently  continuing  his  work  along  such  lines  and  with  some 
promise. 

Upon  rearranging,  equation  (21)  can  be  written  in  the  form* 


2Y|2 

2 

o,.a 


(1.2853)  F  (vi,Gi;  p). 


(28) 


where  F  is  a  function  of  the  material  constants  and  the  angle  P  of  the 
debonded  interface.  A  plot  of  this  equation  for  conditions  of  plane  stress,  as 
well  as  of  plane  strain,  is  given  in  Fig.  8.  In  both  cases  the  maximum  occurs  at 
P  =  60°.  Beyond  this  angle,  the  crack  will  gradually  curve  away  from  the 
interface  and  into  the  matrix. 

In  order  for  us  to  obtain  an  estimate  for  the  critical  stress  for  crack 
initiation  we  let  P  0+,  ie  very  small  but  not  zero.  Thus,  for  our  example 

(Oo)ct  V2a  P  »  1.8186  •^Y|2^1 '  atz  =  0.  (29a) 

On  the  other  hand,  in  the  neighborhood  of  the  free  surface,  the  applied  stress 
is  much  higher  because  of  the  singularity  presence.  In  order  to  overcome  this 
difficulty,  one  may  average  the  local  stress  over  a  distance  equal  to  10%  of  that 
of  the  radius,  ie. 


*  It  should  be  noted  that  at  the  crack  ends  the  stress  field  oscillates  and  that  some  overlap  of 
the  crack  faces  takes  place.  This  matter  is  well  recognized  and  has  been  documented  by 
Williams  (1952),  Rice  et  al.  (1965)  and  England  (1965).  The  region  where  this  occurs,  however, 
is  so  small  (less  than  a  x  lO*^)  that  eq.  (28)  provides  a  good  approximation. 
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(30) 


/  V  1  ro.la  /pn-0.2489 

(4.0633 ao)  [I j 

»  9.5958  Oo  • 

Thus* 

(9.595800)0.  V2aP  -  1.8146  ^7^2 '  atz«±h.  (29b) 


Combining  next  eqs.  (29a)  and  (29b)  one  finds 

(oo)cr  lat  2=h 

-  =  0.10, 

(oo)cr  lat  z=0 


(31) 


ie.  the  critical  loading  stress  which  may  cause  failure  close  to  a  free  surface  is 
approximately  (1/10)  of  the  critical  stress  required  to  cause  the  same  failure  at 
the  center  of  the  fiber's  length.  Thus,  all  things  being  equal,  a  crack  will 
initiate  at  the  free  surface  and  will  propagate  along  the  periphery  of  the 
fiber/matrix  interface  as  well  as  parallel  to  the  axis  of  the  fiber. 

Focusing  next  our  attention  on  the  advancement  of  the  crack  along  the 
periphery  of  the  fiber  we  conclude  that  the  crack  will  advance  itself  to  a 
cricital  angle  of  P  »  60®.  Once  the  crack  has  reached*  8  =60°,  the  local 
geometry  is  similar  to  that  of  a  hole.  This  problem  has  also  been  investigated 
for  the,  3D,  stress  field  close  to  a  ft-ee  surface  (Folias,  1987),  as  well  as  in  the 
interior  of  the  plate  (Folias  and  Wang  1986).  Without  going  into  the  details, 
at  z=h,  it  was  found  that 


suggesting,  therefore,  that  the  failure  now  is  governed  by  the  stress  00^  which 
attains  its  maximum  value  at  8  =  7c/2.  Thus,  the  crack  will  begin  to  curve  into 

*  The  reader  may  r\otice  that  the  right  hand  side  of  equation  (29b)  differs  from  (29a)  because  it 
is  based  on  plane  stress. 
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the  matrix  until  its  direction  becomes  perpendicular  to  that  of  the  applied 
load. 

PERIODIC  ARRAY  OF  FIBERS 

The  previous  results  were  based  on  the  presence  of  one  fiber  only.  It  is 
now  desirable  to  extend  these  results  to  also  include  a  doubly  periodic  array  of 
fibers  which  are  embedded  into  a  matrix.  For  this  reason,  we  assume,  a 
periodic  arrangement  of  the  type  shown  in  fig.  9.  Following  the  same  method 
of  solution  as  that  of  Penado  and  Fblias  (1989),  one  finds*  at  z«0  the  stresses 
<41>and<4e’  ,  for  VI  s  0.34,  V2  «  0.22  and  various  (G2/G1)  ratios,  shown  in  figs. 
6  and  10.  Two  observations  are  worthy  of  note.  First,  beyond  a  certain  ratio  of 
(G2/G1)  the  stress  reaches  an  asymptotic  value.  Such  trend  was  also 

found  by  Adams  (1967)  based  on  2D  considerations.  Second,  as  the  volume  of 
fibers  increases  the  stress  decreases  by  as  much  as  40%  (see  fig.  11). 

Returning  next  to  the  strain  energy  release  rate,  equation  (21)  is  still  a 
good  approximation  provided  that  o©  is  replaced  by  the  following  effective 
load  stress 


f  „a)|^  ]  ■ 

<°o>effeciive  periodic  ’  jiTi  °o  =  W)  oo  , :  for  i  -  0,  (33) 

'X-OJ 

Thus  equation  (29a)  now  becomes 

{,o)„  a  (34) 

which  is  valid  for  small  values  of  B- 


*  The  results  are  valid  for  all  fibers  whidi  are  at  least  four  diameters  away  from  the  bounding 
planes  x  >  ±  w  and  y  »  ±/  The  solution  and  the  details  are  similar  to  those  discussed  by  Penado 

and  Folias  (1989)  except  that  one  now  has  cos{2ne),  n»  0,1  A.-,  where  the  remaining  unknown 
coefficients  are  determined  from  the  boundary  conditions  of  the  geometrical  cell 
configuration.  The  present  results,  are  based  on  n  »  0,...,  N  «  20  terms  which  provide  accurate 

results  in  the  region  I  z/h  I  <  1/2.  However,  many  more  terms  are  needed  in  order  to  obtain 

accurate  results  particularly  in  the  neighborhood  of  z  «  ±  h.  We  are  presently  working  on  this 
and  the  results  for  this  problem,  as  well  as  for  the  problem  of  stresses  due  to  temperature 
missmatch,  will  be  reported  in  the  near  future. 
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Unfortunately,  in  order  to  obtain  a  similar  expression  for  z=h,  one  needs  to 
establish  whether  the  order  of  the  singularity  strength  increases  as  adjacait 
fibers  approach  the  fiber  in  question.  In  view  of  some  previous  work  the 
author  conjectures  that  this  may  very  well  be  the  case.  Thus,  the  following 
fundamental  questions  come  to  mind.  How  close  must  adjacent  fibers  be 
before  the  order  of  the  singularity  strength  is  affected?  Does  a  certain 
separation  distance  or  a  certain  periodic  array  of  fibers  exists  which  leads  to  an 
optimal  state  of  stress?  Based  on  3D  considerations,  Penado  and  Folias  (1989) 
have  shown  that  when  fibers  are  placed  four  fiber  diameters  apart,  center  to 
center,  practically  all  fiber  interactions  have  subsided,  including  those  at  the 
free  surface  2=h.  The  author  suspects,  however,  that  when  fibers  are  placed 
two  diameters  apart,  center  to  center,  the  singularity  strength  will  be  affected. 
Natxirally,  this  is  a  conjecture  that  needs  to  be  investigated. 

As  a  practical  matter,  if  one  uses  the  approximation  given  by  eq. 

(30),  the  critical  stress  to  failure  at  the  fiber  edge  in  a  glass 
fiber/epoxy  matrix  composite  with  the  properties 

Gm  —  2.lOGPa  Vn  0.34  a/  ss  10"®cm  =  60® 

G/*35.00GPa  i//a0.22  27„  70J/m’  V/ *  0.70 

becomes 

(oo)cr  —  20.582  F(Vf)  Mpa. 

=  2.985  F(Vf)  ksi  ;  at  the  fiber  edge 

a  plot  of  which  is  given  in  Fig.  12.  Edge  delamination  may  now  be 
modeled  as  the  progressive  failure  of  a  row  of  fibers. 


(35) 


(36) 
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CONCLUSIONS 


Based  on  a  3D  analytical  solution,  we  have  shown  that  fiber/matrix 
debonding  is  most  likely  to  occur  dose  to  a  free  surface.  Thus,  regions  where 
fibers  intersect  free  surfaces,  eg.  holes,  cut  outs,  edges,  cracks  etc.  are  potential 
trouble  spots.  Moreover,  the  strain  energy  release  rate  (eq.  28)  may  be  used  to 
predict  crack  initiation  in  the  center  of  fl\e  fiber  length  (eq.  29a),  as  well  as  at 
the  free  surface  (eq.  29b).  Moreover,  fiber/matrix  debonding  at  a  free  surface 
will  occur  at  approximately  (1/10)  the  load  value  required  for  the  same  type  of 
failure  to  occur  at  the  center  of  the  fiber  length.  Such  information  on  crack 
initiation  is  particularly  important  for  the  proper  understanding  of  damage 
evolution. 

Alternatively,  the  strain  energy  release  rate  for  a  periodic  array  of  fibers 
of  the  t3rpe  shown  in  fig.  9  may,  at  2*0,  be  approximated  by  eq.  (28)  in 
conjunction  with  eq.  (33).  A  similar  expression  applicable  to  the 
neighborhood  of  the  free  surface  requires  that  one  must  first  establish 
whether  the  strength  of  the  singularity  is  indeed  affected  as  the  fiber  volume 
increases.  For  Vf  ^  0.05,  however,  it  has  been  shown*  that  no  such  interaction 
effects  are  present. 

As  a  final  remark,  we  note  that  if  the  bond  at  the  interface  does  not  fail 
the  analysis  shows  that  there  exists  a  stress  magnification  factor  in  the  resin 
which  attains  a  maximum  between  the  fibers.  This  maximum  stress 
magnification  occurs  along  the  line  6  »  0^  and  at  a  distance  r  *  1.2a  from  the 
center  of  the  fiber** . 


*  See  Panado  and  Folias  (1989). 

**  This  condition  is  valid  for  all  0  S  z  <  h. 
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HGURE  CAPTIONS 


Fig.  1.  Geometrical  and  loading  configuration. 

Fig.  2.  Definition  of  local  coordinates. 

Fig.  3.  Singularity  strength  for  isotropic  fiber  and  isotropic  matrix 
versus  G2/G1. 

Fig.  4.  Stress  Ojy^at  r*a,  0=0  and  for  vi  =  0.34,  V2  =  0.22  and  (G2/G1)  = 
16.67,  accross  the  thickness. 

Fig.  5  Stress  at  r=a,  0=0  and  for  vi  =  0.34,  V2  =  0.22  and  (G2/G1)  = 

16.67,  across  the  thickness. 

Fig.  6  Stress  oJJ^at  r=a,  0=0  and  for  vi  =  0.34,  V2  =  0.22,  versus  the  ratio 

(G2/G1). 

Fig.  7  Fiber /matrix  interface  crack  under  transverse  loading. 

Fig.  8  Strain  energy  release  rate  for  plane  stress  and  plane  strain 

conditions  for  vi  =  0.34,  V2  =  0.22  and  (G2/G1)  =  16.67. 

Fig.  9  Periodic  array  of  fibers  of  length  2h,  embedded  into  a  matrix. 

Fig.  10  Stress  oJJ^at  r=a,  0=0  and  for  vi  =  0.34,  V2  =  0.22,  versus  the  ratio 

(G2/G1). 

Fig.  11  Stress  at  r=a,  versus  Vf,  for  (G2/G1)  =  16.67  vi  =  0.34,  V2  = 

0.22. 

Fig.  12.  Critical  stress  versus  Vf  for  (G2/G1)  =  16.67,  v]  =  0.34  and  V2  = 

0.22. 
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Fig.  2. 


Definition  of  local  coordinates. 
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Fiber/matrix  interface  crack  under  transverse  loading. 


STRAIN  ENERGY  AND  ITS  RELEASE  RATE 
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Abstract 


In  this  paper,  the  3D  stress  field  is  derived  for  a  cylindrical  fiber  consisting  of  a 
radius  a  and  surrounded  by  a  itKxlified  cylindrical  shell  matrix  of  outside  radius  b.  The 
system  is  embedded  into  a  plate  matrix  of  arbitrary  thickness  2h.  All  materials  are  assumed 
to  be  isotropic  and  perfect  bonding  is  assumed  to  prevail  at  their  respective  interfaces.  As 
to  loading,  a  uniform  tension  is  applied  in  the  plane  of  the  matrix  plate  and  at  points  remote 
from  the  interfaces. 

The  analysis  shows  that  all  stresses  are  sensitive  to  the  radius  to  half  thickness  ratio 
(aA))  as  well  as  the  material  properties.  In  general,  the  effect  that  a  modified  shell  matrix 
has  on  the  magnitude  of  the  radial  stress  Ojr  is  minimal  (approximately  3%).  However,  the 
location  of  its  max.  depends  strongly  on  the  respective  shear  moduli  ratios  (G2A3i)  and 
(G3/G1).  The  analysis  also  characterizes  the  type  of  fracture  process  to  be  primarily 
cohesive.  Moreover,  the  stresses  are  shown  to  be  singular  in  the  neighborhood  of  the 
intersection  of  the  interface  curve  with  that  of  the  edge  plane. 

All  things  being  equal,  a  crack  is  most  likely  to  originate  at  the  free  edge  and  will 
propagate  ,  along  the  interface  and  towards  the  interior  of  the  composite  ,up  to 
approximately  one  fiber  radius.  Whereby  it  will  ,subsequently,  curve  into  the  matrix  and 
follow  the  locus  of  the  max.  On.  Approximate  failure  criteria  arc  also  suggested. 


INTRODUCTION 


The  properties  of  fiber  reinforced  composites  are  very  much  dependent  upon  the 
stability  of  the  interfacial  region  between  the  matrix  and  the  fiber .  This  is  because  the 
primary  function  of  this  interface  region  is  to  transmit  a  portion  of  the  load  from  the  matrix 
to  the  reinforcing  fibers  and  vise  versa.  This  ability  to  transmit  stress  across  the  phase 
boundary  must ,  ultimately ,  depend  upon  the  mechanical  properties  of  the  matrix,  the  load 
bearing  capacity  of  the  fibers  as  well  as  the  strength  of  the  fiber  matrix  interface.  It  is 
natural  .therefore,  to  seek  the  relationship  between  the  overall  composite  strength  and  the 
above  variables. 

In  a  recent  paper  which  was  based  on  3D  considerations,  Panado  and  Folias  (1989) 
studied  the  problem  of  a  cylindrical  (isotropic)  fiber  which  is  embedded  into  an  epoxy 
matrix.  The  interface  between  the  fiber  and  the  matrix  was  represented  as  a  cylindrical 
boundary.  In  reality,  however,  it  should  be  an  interface  region  rather  than  an  interface 
boundary  the  composition  of  which  depends  on  the  specific  fiber/matrix  combination 
(Drzal,  1983).  There  are  three  reasons  for  this.  First,  there  exist  morphological  variations 
on  the  surface  of  the  fiber.  Second,  there  is  a  chemical  reaction  which  takes  place  in  such 
regions.  Third,  in  order  to  protect  the  fibers  from  mechanical  damage  and  to  promote  high 
matrix  adhesion,  mixtures  of  sizing  and  coupling  agents  are  applied  to  the  fiber  surfaces 
during  the  fiber  manufacturing  process.  The  question  ,  therefore,  arises  as  to  what  effect  if 
any  does  such  an  imperfect  interface ,  or  modified  matrix  layer ,  have  on  the  strength 
characteristics  of  the  overall  system.  In  oracr  to  provide  a  partial  answer  to  this,  Agarwal 
and  Bansal  (1979)  investigated  the  effects  of  such  an  'imperfect  interface’  by  considering  a 
thin  layer  of  different  material  properties  to  represent  the  interface  region.  The  analysis  wa.s 
restricted  to  a  load  in  the  longitudinal  direction  of  the  fiber  and  was  based  on  2D 
considerations.  Subsequently,  Tandon  and  Pagano  (1988)  developed  an  approximate 
model  that  one  may  use  to  predict  the  overall  moduli  and  coefficients  of  thermal  expansion 


of  composites  in  the  presence  of  imperfect  bonding  and  subject  to  displacement  and  traction 
boundary  conditions. 

The  present  paper  investigates  the  effect  that  such  a  'modified  layer’  has  on  the 
strength  characteristics  of  an  isotopic  fiber  embedded  into  a  matrix  and  subjected  to  a 
uniform  transverse  load.  Moreover,  in  order  to  extract  the  relevant  edge  effects  the 
analysis  will  be  based  on  3D  considerations .  It  is  hoped  that  this  study  will  compliment 
the  previous  investigations  and  perhaps  offer  some  further  insight  on  the  possibilities  of 
cracking  due  to  the  external  load. 

FORMULATION  OF  THE  PROBLEM 

Consider  the  equilibrium  of  a  body  which  occupies  the  space  Ixl  < «»,  lyl  <  oo, 

Izl  <  h  and  contains  three  regions  with  different  elastic  propenies.  Their  common 
boundaries  consist  of  through-the-thickness  cylindrical  surfaces  of  radii  r  =  a  and  r  =  b, 
whose  generators  are  parallel  to  the  2-  axis  (see  Figs  1  and  2).  The  situation  described  is 
that  of  a  fiber,  denoted  by  superscript  (3),  which  is  embedded  into  a  matrix,  denoted  by 
superscript  (1).  Moreover,  the  fiber  surface  is  coated  with  a  thin  layer  of  different  material 
propenies,  denoted  by  superscript  (2).  All  three  materials  are  assumed  to  be 
homogeneous,  isotropic  and  linear  elastic.  At  the  interfaces  r  =  a  and  r  =  b  perfect  bonding 
is  assumed  to  prevail.  As  to  loading,  a  uniform  tensile  stress  Oq  is  applied  at  the  matrix  at 
points  remote  from  the  fiber.  In  all  three  regions,  the  surfaces  Izl  =  h  are  assumed  to  be  free 
of  stress  and  contstrain  . 


In  the  absence  of  body  forces,  the  coupled  differential  equations  governing  the 
displacement  functions  u  are; 


JUl 

mi-2 


3eW 


+  V^ 


i,  k  =  1, 2,  3 


(1) 


'  ‘jj  «ik  ^ '  'ft’}  2. 3 

where  are  the  respective  shear  moduli. 

As  to  boundary  conditions,  we  require  that: 


(2) 


(3) 


as  Ixl  -» <»  :  o  =  1  =  1  =  0 

as  lyl  :  t  =  i  =  0,  o  =  o„ 
at  Izl  =  h :  t^i,’  =  X  ®  =  o  ®  =  0  ;  i  =  1, 2, 3 

a.r  =  b:  o<i)  -  o®  =  X®-  X®  =  x®-  x®  .0, 

rr  rr  rw  ru  rz  rz 
at  r  =  b  ;  u  -  u  =  v  -  v  =  w  -  w  =  0  , 

atr  =  a-  (3)_  (2)_  ^  (3)  ^  (2)  _  (3) 

atr-a.  ^TT~^r0  ^rO  ^rz  Sz"^’ 

at  r  =  a  ;  u  -  u  =  v  -  v  =  w  -  w  =  0  , 

r  r  o  o 


Finally,  at  r  =  0  we  must  require  that  all  stresses  and  displacement  be  bounded. 


(4) 

(5) 

(6) 

(7) 

(8) 
(9) 

(10) 


METHOD  OF  SOLUTION 

A  general  soludon  to  a  certain  class  of  three  dimensional  boundary  value  problem.s 
which  arise  in  elastostatics,  was  developed  by  Folias  (1976)  and  was  subsequently  put  in  a 
much  more  convenient  form  (Folias  1990).  The  general  solution  was  recently  used  to 
solve  for  the  3D  stress  field  of  one  cylindrical  fiber  embedded  into  an  infinite  matrix 


(Penado  and  Folias  1989).  The  latter  study  may  now  be  extended  to  also  account  for  the 
effects  of  a  modified  shell  matrix  surrounding  the  fiber. 

The  mathematical  details,  however,  are  long  and  tedious  but  similar  to  those 
of  the  above  reference,  except  that  one  now  has  to  deal  with  two  interfaces  instead  of  one. 
Thus  for  the  sake  of  brevity,  the  details  in  the  present  analysis  are  omitted  .  In  summary, 
the  displacement  fields  are  expressed  in  terms  of  two  infinite  series,  one  with  complex 
eigenfunctions  and  the  other  with  real  eigenfunctions.  By  construction,  the  general 
solution  automatically  satisfies  the  free  of  stress  boundary  conditions  at  Izl  =  h.  The 
unknown  series  coefficients  are  subsequently  determined,  numerically,  from  the  remaining 
boundary  conditions  at  the  two  interfaces  r  =  a  and  r  =  b.  For  more  details,  the  reader  is 
referred  to  the  work  by  Penado  and  Folias  (1989). 

NUMERICAL  RESULTS 

Once  the  series  coefficients  have  been  determined,  the  3D  displacement  and  stress 
fields  may  then  be  recovered.  Throughout  the  numerical  work,  we  have  used  double 
precision  because  the  system  is  very  sensitive  to  small  changes.  Furthermore,  the  rate  of 
convergence  has  been  addressed  in  the  work  by  Penado  and  Folias  (1989).  The  present 
numerical  work  also  exhibits  similar  convergence  characteristics.  Thus,  without  going 
into  the  numerical  details,  the  displacement  w ,  on  the  surface  z  =  h  ,  as  a  function  of 

is  shown  in  Figs.  3  and  4  for  various  G  ratios.  This  particular  displacement  is  of  special 
interest  to  the  experimentalists  for  future  comparisons. 

The  behavior  of  the  radial  stress  Orr,  as  a  function  of  the  modified  layer  thickness 
A  =  (b-a)/a,  is  shown  in  Figs  5  and  6  ,  at  the  position  0  =  0  and  r  =  b.  It  is  noted  that  for 
the  fixed  ratios  of  (G3/G1)  =  16.67  and  (G2/G1)  =  12,  the  relative  max.  of  the  radial  stress 
occurs  at  A  =  0.2  ,  for  z  =  0  and  z  =  h.  In  fact,  the  same  trend  prevails  throughout  the 
thickness  z.  In  general,  the  location  of  the  max.  is  a  function  of  the  respective  material 


properties.  Moreover,  the  behavior  of  the  radial  stress  On-  versus  (z/h)  is  similar  to  that 
found  by  Penado  and  Folias  (1989).  A  typical  plot  is  shown  in  Fig.  7.  On  the  other  hand, 
the  variation  of  the  stress  Oee  as  a  function  of  A ,  at  6  =  0  and  r  =  b  ,  is  given  in  Figs.  8 
and  9  for  z  =  0  and  z  =  h  .respectively.  Thus  as  A  increases,  Oqq  decreases  slightly  (by 
approximately  3%). 

Fig.  10,  depicts  the  variation  of  the  radial  stress  0^  as  a  function  of  the  ratio 
(G3/G2)  for  6  =  0,  z  =  0,  A  =  0.2,  r  =  b  and  (G3/G1)  =  16.67.  Similarly,  Figs.  1 1  and  12 
show  the  variation  of  On- .  for  (G2/G1)  <1  and  for  (G2/G1)  >  1  ratios,  respectively.  Here 
again  the  max.  change  is  noted  to  be  approximately  3%. 

It  is  well  recognized  from  experimental  observations  that  debonding  along  a 
fiber/matrix  interface  is  a  common  occurrence  in  laminated  composite  plates  due  to 
transverse  loads.  What  is  not  clear,  however,  is  whether  the  crack  originates  at  the 
fiber/matrix  interface  or  in  the  matrix  itself  but  close  to  the  interface.  Naturally,  all  things 
being  equal,  a  crack  is  most  likely  to  originate  at  a  location  where  the  radial  stress  attains  its 
maximum.  To  provide  us  with  some  further  insight  on  this  matter,  we  plot  in  Fig.  13  the 
radial  stress  On-  versus  the  distance  (r-a) ,  for  various  (G2/G1)  ratios.  The  reader  may 
notice  that,  although  the  stress  On-  is  continuous  at  r  =  b,  its  derivative  there  is 
discontinuous.  This  result  meets  our  expectations  for  the  material  properties  to  the  left  and 
to  the  right  are  distinct.  Furthermore,  it  may  be  noted  that  for  (G2/Gi)<  2.5  ratios  the 
max.  radial  stress  occurs  within  the  modified  shell  matrix,  while  for  2.5  <  (G2/G1)  occurs 
within  the  matrix,  i.e.  past  the  interface  r  =  b.  Thus  the  question  as  to  whether  one  has 
adhesive  or  cohesive  fracture  really  depends  on  the  respective  material  properties.*  It 
should  be  emphasized  that  this  comment  is  applicable  not  only  to  the  case  of  one  fiber ,  but 
also  to  the  case  of  an  infinite  periodic  array  of  fibers  which  are  spaced  by  four  (or  more  ) 
fiber  diameters  center  to  center.  In  the  laner  case,  there  exist  no  significant  interactions 
present  between  the  adjacent  fibers  (see  Figs.  3  and  4).  This  matter  has  also  been 

♦From  Fig.  13,  it  is  noted  that  away  from  the  edge,  adhesive  fracture  only  takes  place  at 
one  value  which  is  between  1.7  <  (G2/G1)  <2.8. 


discussed  in  a  recent  paper  by  Folias  (1991)  which  is  based  on  3D  considerations.  In  the 
event  that  the  fibers  are  closer  together,  Adams  et  .al .  (  1967  )  have  shown,  on  the  bases 
of  2D  considerations,  that  the  max.  radial  stress  0^  occurs  at  r  =  1.2a,  i.e.  between  the 
fibers.  The  same  result  has  subsequently  been  confirmed  by  other  researchers  in  the  field. 
On  the  bases  of  3D  considerations  and  the  given  data ,  the  authors  also  confirm  this  result 
(Liu  and  Folias  1991)  which  is  valid  throughout  the  interior.  In  general,  given  a  set  of 
material  properties,  the  location  of  the  maximum  is  a  function  of  the  ratio  (a/b).  Returning 
next  to  the  free  edges,  the  stress  field  in  the  neighborhood  of  the  intersection  of  the 
interface  boundary  with  that  of  the  free  surface  was  shown  to  be  singular  (Folias  1989,  Li 
and  Folias  1991).  The  same  result  is  expected  to  apply  here  too.  In  fact,  the  authors 
believe  that  the  strength  of  the  stress  singularity  will  increase  slightly  as  the  distance 
between  the  adjacent  fibers  decreases  and  the  interaction  effects  become  stronger.  Similar 
trends  are  also  expected  to  prevail  in  the  case  of  the  modified  layer  matrix. 

CONCLUSIONS 

The  following  conclusions  may  be  drawn  from  the  previous  results: 

(i)  away  from  the  free  edge: 

•  the  effect  of  a  modified  shell  matrix  on  the  magnitude  of  the  radial  stress  Orr  is 
negligible,  ie.  approximately  3% 

•  the  max.  value  of  the  radial  stress  On  is  approximately  the  same  for  all 
(G2/G1)  ratios  examined 

•  the  location  of  this  max.  occurs  within  the  modified  layer  if  1  <  (G2/G1)  <  2  , 
and  within  the  matrix  if  2  <  (G2/G1)  <  16.67 

•  characterization  of  the  fracture  process  as  adhesive  or  cohesive  depends  on  the 
material  properties  (see  Fig.  13) 

•  the  modified  layer  thickness  A  which  will  make  the  radial  stress  On  a  relative 
maximum  is  .20 

(ii)  within  the  vicinity  of  the  free  edge: 

•  the  max.  radial  stress  On  occurs  at  one  of  the  two  interfaces 


In  conclusion,  all  things  being  equal,  a  crack  is  most  likely  to  originate  in  the 
vicinity  of  the  free  edge  and  along  one  of  the  interface  boundaries  (that  of  the  highest  stress 
singularity).  Subsequently,  the  crack  will  propagate  away  from  the  edge  and  along  the 
interface  up  to  a  distance  of  one  fiber  radius  at  which  point  it  will  slowly  curve  into  the 
manix  and  follow  the  locus  of  the  max.  radial  stress  On-  (see  Fig.  13).  Thus,  it  is  now 
possible  to  derive  two  separate  failure  criteria  which  will  be  applicable  to  the  regions  (i) 
close  to  the  free  edge  and  (ii)  away  from  the  free  edge.  The  first  author  has  already  derived 
such  a  criterion  for  the  prediction  of  failure  at  a  free  edge  (see  Folias,  1991) .  It  remains, 
therefore  ,  for  us  to  derive  a  criterion  applicable  to  the  latter  region. 
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Geometrical  configuration. 


Fig.  2 


Geometrical  omfiguration  at  the  iree  edge. 
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The  displacement  w  at  the  free  surface  z  =  h  vs.  (r/a)  for  (a/h)  =  4. 
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Fig.  6  The  radial  stress  vs.  A  at  r^b  and  h. 
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The  max.  circumfcrcmial  stress  vs.  A  at  z=  h. 


Fig.  1 1  The  max.  radial  stress  vs  the  ratio  (r-a)/a  for  various  G2/G1  ratios  at  z=0. 
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ABSTRACT 


A  three-dimensional  (3-D)  micromechanical  model  has  been  developed  to  represent  a 
unidirectional  composite  plate  that  is  subjected  to  a  residual  stress  load  induced  due  to  a 
temperature  mismatch.  The  model  assumes  the  fibers  to  be  cylindrical  inclusions  that  are 
periodically  embedded  into  an  epoxy  matrix  plate.  At  the  fiber/matrix  interface,  perfect 
bonding  conditions  are  assumed  to  prevail.  The  method  used  in  the  thesis  is  based  on  a  3-D 
analytic  considerations. 

The  analytical  solution  shows  the  maximum  stresses  to  occur  at  0=0  for  the  case  of  a 
graphite/epoxy  composite  plate.  The  location  of  the  maximum  of  the  stresses,  however, 
may  change  for  other  material  systems,  particularly,  for  small  shear  moduli  ratios 
The  solution  also  shows  that  the  shear  moduli  ratio  \3i2j\i\  (fiber  to  matrix),  the  geometric 
ratios  a/h  (fiber  radius  to  thickness  of  the  plate)  and  a/b  (fiber  radius  to  separation  distance 
between  adjacent  fibers)  effect  greately  the  displacement  and  stress  distributions.  On  the 
other  hand,  the  ratio  of  thermal  expansion  coefficients  a(2)/a(i)  (fiber  to  matrix)  presents  a 
linear  action  to  the  displacement  and  stress  fields.  The  stress  profiles  along  a  fiber  length 
are  shown  to  be  constant  except  in  the  neighborhood  of  the  fiber  edge  where  a  boundary 
layer  is  shown  to  prevail.  In  this  region,  a  stress  singularity  is  shown  to  prevail  which  is  a 
departure  from  the  results  obtained  by  macromechanical  theories. 

In  the  limit,  as  a/b-»0,  the  result  of  one  fiber  model  is  recovered;  as  a/h— >oo,  the  plane 
stress  solution  is  recovered.  Similarly,  as  the  result  of  a  plate  with  a  periodic 

array  of  holes  is  recovered. 


CHAPTER  1 


INTRODUCTION 

It  is  well  recognized  that  fiber-reinforced  composite  materials  are  very  attractive  for 
their  high  strength  and  stiffness  to  weight  ratios  and  for  their  excellent  fatigue  resistance 
characteristics.  With  the  advancement  of  composites  and  their  extensive  use  in  aerospace, 
automotive  and  other  industries,  problems  for  the  determination  of  the  displacement  and 
stress  fields  have  drawn  considerable  attention.  However,  the  anisotropic  behavior  of 
composite  materials  makes  the  stress  field  more  complex  than  that  of  traditional  engineering 
materials.  Many  designers  encounter  difficulties  in  determining  the  stress  field  around  the 
interface  of  fiber  and  matrix  and  a  free  surface.  Thus  reliable  prediction  of  the  stress  field 
within  a  composite  is  an  important  research  topic  which  needs  further  investigaticMi. 

Basically,  there  are  two  categories  of  factors  that  may  induce  materials  failure,  the  first 
one  is  due  to  the  external  loads  and  the  second  is  due  to  localized  stress  fields,  for  example 
residual  stresses  caused  by  the  environment  such  as  a  temperature  change.  Upon  the 
application  of  external  loads,  regions  of  low  toughness  values  may  fail  before  the  bulk  of 
the  material  does  due  to  the  presence  of  impurities,  inclusions,  grain  boundaries,  etc.  On 
the  other  hand,  the  residual  stresses  are  present  as  a  result  of  the  manufacturing  process  of 
the  material.  During  the  change  of  temperature,  residual  stresses  are  induced  because  of  a 
thermal  expansion  mismatch,  thermal  expansion  anisotropy,  or  phase  transformation. 
Depending  on  the  values  of  the  thermal  expansion  coefficients,  the  residual  stresses  may  be 
tensile  or  compressive  in  nature  and  both  types  of  stresses  may  induce  microcracking. 

In  order  to  be  able  to  predict  the  failing  characteristics  of  composite  materials, 
particularly  in  the  neighborhood  of  a  free  surface  and  particularly  at  the  fiber/matrix 
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interface,  it  is  necessary  to  know  the  local  stress  field  fixmi  a  3-D  point  of  view.  Numerous 
researchers  have  focused  on  the  effects  which  one  inclusion  has  in  a  composite.  In  1942, 
Reissner  [1]  introduced  a  semi-direct  variational  method  consisting  of  3-D  corrections  to 
the  theory  of  plane  stress.  In  1948,  Green  (2]  proposed  a  3-D  solution  in  terms  of  infinite 
series  of  a  very  complicated  nature.  His  solution  satisfied  the  governing  differential 
equations,  but  the  boundary  conditions  were  satisfied  by  an  approximate  iterative 
procedure.  One  year  later.  Green  [31  developed  a  general  method  based  on  infinite  series 
for  the  solution  of  3-D  boundary  value  problems.  However,  the  form  of  his  solution  is 
rather  complicated  for  practical  use.  Alblas  [4],  in  1957,  used  the  general  procedure 
proposed  by  Green  to  construct  an  infinite  series  representation  of  the  solution  ’n  a  simple 
form.  He  expressed  the  solution  to  Navier's  equations  in  terms  of  a  set  of  complex 
eigenfunctions.  In  1961,  Reiss  [5]  sought  the  correction  factor  to  the  plane  stress  solution 
on  the  basis  of  an  asymptotic  expansion  of  the  stress  field  in  powers  of  the  thickness 
parameter.  Two  years  later,  he  applied  his  theory  to  find  a  3-D  correction  to  Kirch’s 
solution.  The  results,  however,  are  valid  only  for  relatively  thin  plates.  Youngdahl  et. 
al.[6],  in  1966,  presented  an  integral  form  of  the  solution  by  means  of  a  specially  adapted 
integral  transform  technique.  Their  solution,  however,  is  only  applicable  to  the  case  of  a 
half  space.  Later,  in  1975,  Folias  [7]  developed  a  method  for  solving  3D  mixed-boundary- 
value  problems  that  arise  in  elastostatics.  The  method  is  applied  to  a  plate  of  finite 
thickness,  which  contains  a  finite,  through  the  thickness,  line  crack.  Based  on  Folias 
previous  woik,  Folias  and  Wang  (1990)  specialized  the  general  solution  to  the  equilibrium 
of  a  linear  elastic  layer  for  the  case  of  a  plate  of  an  arbitrary  thickness  that  has  been 
weakened  by  a  cylindrical  hole  [8].  The  solution  satisfied  both  Navier’s  equations  as  well 
as  the  plate  boundary  conditions.  Moreover,  in  1989  Folias  [9]  investigated  analytically  the 
stress  field  in  the  neighborhood  of  the  intersection  of  a  hole  and  a  free  surface.  The  analysis 
shows  the  complementary  solution  to  be  proportional  to  p®,  where  p  is  the  distance  ftom 
the  point  of  intersection  of  the  hole  and  the  free  surface  of  the  plate.  The  analysis  showed 
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that  no  stress  singularity  was  present  in  this  vicinity.  In  1989,  Penado  and  Folias  [10] 
extended  the  previous  work  to  the  case  of  a  cylindrical  inclusion.  Both  the  plate  and  the 
inclusion  are  assumed  to  be  of  homogeneous  and  isotropic  materials  with  different  material 
properties.  In  the  limit,  as  the  ratio  of  shear  moduli  =  0.00001,  the  results  for  a 
cylindrical  hole  are  recovered.  At  the  same  time,  Folias  [11]  studied  analytically  the  stress 
fleld  in  the  neighborhood  of  the  intersection  of  the  cylindrical  inclusion  and  the  free 
surface.  The  displacement  and  stress  fields  were  derived  explicitly  and  the  presence  of  a 
stress  singularity  was  shown  to  exist.  Moreover,  this  stress  singularity  is  shown  to  be  a 
function  of  the  ratios  of  the  shear  moduli  and  Poisson’s  ratios.  In  1990,  Folias  and  Liu 
[12]  extended  Penado  and  Folias  '  work  to  a  more  general  case  which  includes  a  coating 
layer  between  the  fiber  and  the  matrix.  In  1991,  Zhong  and  Folias  [13]  investigated  the 
loading  transfer  characteristics  between  a  fiber  and  a  matrix  plate  subjected  to  an  axial 
loading.  Both  perfect  bonding  and  imperfect  bonding  conditions  were  examined. 

As  for  periodical  systems,  Goree  [14]  in  1967  presented  a  solution  for  the 
displacements  and  stresses  in  an  infinite  elastic  matrix  containing  two  perfectly  bounded 
rigid  circular  cylindrical  inclusions  of  different  radii,  and  of  infinite  length  normal  to  the  x-y 
plane.  In  the  same  year,  Hedgepeth  [15]  obtained  a  solution  for  two  stress  distribution 
problems  that  resulted  from  breaking  of  the  filaments  in  a  composite  material  composed  of 
high  modulus  elements  embedded  in  a  low  modulus  matrix.  Subsequently,  Adams  (1967) 
presented  two  papers  where  the  longitudinal  shear  loading  [16]  and  transverse  normal 
loading  [17]  were  discussed  respectively.  Using  an  analysis  based  on  the  theory  of 
elasticity  ,  the  problem  of  a  doubly  periodic  rectangular  array  of  elastic  filaments  in  an 
elastic  matrix  material  was  formulated  in  his  papers  where  the  effect  of  a  uniform 
temperature  was  also  imposed.  Later,  Haener  [18]  used  the  displacement  potential  method 
to  express  the  3-D  stress  distributions  in  a  unidirectional  multi-fiber  composite  under 
external  and  residual  loads.  He  assumed  the  geometric  arrangement  of  the  fibers  was 
hexagonal  and  the  hexagon  boundary  remained  regular  under  the  load.  In  1969,  Marloff 
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and  Daniel  [19]  used  a  experimental  method  to  determine  the  3-D  stress  distribution  in  a 
unidirectional  composite  plate  subjected  to  matrix  shrinkage  and  nornud  transverse  load.  In 
1973,  Keer  et  al.  [20]  used  finite  integral  transforms  to  study  the  phenomena  of  the 
separation  of  a  smooth  circular  inclusion  fiom  a  matrix  that  is  subjected  to  a  uniform  load. 
A  2-D  solution  for  plates  with  perfectly  bonded  circular  inclusions  can  be  found  in  the 
papers  [21]  of  Sendeckyj  (1970).  Based  on  Sendeckyj's  previous  work,  Yu  and  Sendeckyj 
(1974)  extended  their  study  to  the  case  of  an  infinite  elastic  matrix  containing  multiple 
inclusions  by  using  the  Schwarz  alternating  method  [22].  Their  results  were  subsequendy 
specialized  to  cases  of  two  and  three  inclusions.  In  1984,  Adams  [23]  used  a  finite  element 
method  to  analyze  a  microscopic  region  of  a  unidirectional  composite  plate  in  plane  strain 
case.  In  1985,  Evans  et  al.  [24]  analyzed  matrix  fracture  in  britde-matrix  fiber  composites. 
The  numerical  solution  shows  that  matrix  cracks  initiate  prior  to  fiber  failure  and  the 
influence  of  the  fibers  that  bridge  the  matrix  crack  is  represented  at  the  crack  surfaces.  In 
1990,  based  on  2-D  consideration,  Isida  [25]  used  a  complex  series  approach  to  solve  the 
problem  of  an  array  of  periodic  inclusions  embedded  into  an  infinite  matrix  subjected  to  a 
uniform  transverse  load.  Thus,  his  contribution  provides  us  with  further  insight  into  the 
strength  of  the  composite  plate.  In  1991,  Folias  and  Liu  extended  the  previous  works  from 
one  fiber  to  an  array  of  periodic  fibers.  Two  investigations  have  been  completed  for  the 
model  of  a  matrix  plate  embedded  with  a  periodic  array  of  fibers  subject  to  different 
loadings.  One  is  subjected  to  a  uniform  transverse  loading  [26],  the  other  is  subjected  to  a 
axial  loadings  [27].  The  results  show  that  the  effect  of  the  fiber  volume  fraction  Vf  is  an 
important  parameter  that  affects  greatly  the  displacement  and  stress  fields.  In  the  limit,  as 
the  ratio  a/b-^  ( a/b  is  the  ratio  of  the  fiber  radius  to  separation  distance  between  adjacent 
fibers),  the  results  of  Penado  and  the  results  of  Zhong  are  recovered,  respectively. 

In  the  area  of  residual  stress  problems,  however,  most  of  the  research  is  limited  to  a 
macro-approach  which  accounts  only  for  the  macromechanical  behavior  of  the  laminates 
and  neglects  the  individual  behavior  of  the  fibers  and  of  the  matrix.  Griffin  [28]  in  1983 
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used  a  fully  3-D  finite  element  analysis  to  predict  thermal  stress  distributions  in  thick 
graphite/epoxy  laminates.  The  result  shows  that  high  interlaminar  thermal  tensile  stresses 
ate  present  near  the  free  edges.  Later,  in  another  paper  [29],  provided  experimental 
evidence  with  which  he  confirmed  his  numerical  solution.  Another  practical  model  for 
matrix  microcracking  present  in  graphite/epoxy  composites  was  examined  by  Bowles 
(1984)  [30],  where  finite  element  analysis  was  used  to  predict  the  effect  of  different 
thermal  expansion  coefficients.  In  1987,  Evans  et  al.[31]  investigated  the  high-temperature 
mechanical  properties  for  a  ceramic  matrix  crxnposite.  The  change  in  mechanical  behavior 
has  been  attributed  to  a  large  variation  in  the  shear  resistance  of  the  fiber/matrix  interface.  In 
1988,  Garg  [32]  studied  the  fracture  behavior  of  graphite/epoxy  composite  effected  by  the 
change  of  temperature  with  the  experimental  method.  It  proves  that  a  temperature  change 
may  induce  failure  in  the  material.  Fang  [33]  in  1989  presented  an  experimental  and 
analytical  investigation  of  thermally-induced  cracking  in  graphite/epoxy  composites,  where 
the  thermal  stresses  will  affect  the  extent  and  form  of  damage.  However,  his  analysis  is 
based  only  on  an  approximate,  2-D  micro-cracking  model,  which  can  not  be  used  to  predict 
the  stress  distribution  along  the  thickness.  In  1985,  Mikata  et  al.  [34]  calculated  a  stress 
field  in  a  coated  continuous  fiber  composite  subjected  to  thermomechanical  loadings.  The 
investigation  shows  that  the  maximum  stress  occurs  in  the  coating  and  it  reduces  as  the 
volume  fraction  of  fibers  and  coating  thickness  increases.  A  prediction  of  thermal 
expansion  for  unidirectional  composites  was  investigated  by  Bowles  et  al.  [35]  in  1988, 
where  several  analyses  were  compared  with  each  other.  Also  a  sensitivity  analysis  was 
conducted  to  determine  the  relative  influence  of  constituent  properties.  In  1988,  an 
investigation  was  carried  out  by  Delale  [36]  where  the  microcracking  of  ceramic-matrix 
composite  due  to  residual  stresses  was  studied.  Starting  with  a  model  of  single  fiber 
embedded  in  an  infinite  matrix,  the  critical  fiber  size  was  determined  and  the  relationship 
between  critical  fiber  size  and  the  parameters  such  as  thermal  expansion  coefficients  of  fiber 
and  matrix,  change  of  temperature  was  established.  It  is  shown  that  a  small  difference  of 
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the  thermal  expansion  coefficient  leads  to  a  large  critical  fiber  size.  In  1989,  Evans  [37] 
investigated  the  debonding  properties  of  brittle-matrix  composite  subjected  to  residual 
stresses.  The  debonding  behavior  is  shown  to  depend  sensitively  on  the  thermal  expansion 
mismatch.  Later,  in  1990,  Evans  [38]  focused  his  research  on  the  effect  of  the  interface  in 
fiber-reinforced  ceramics.  It  is  shown  that  the  residual  stress  arisen  from  thermal 
contraction  mismatch  upon  cooling  has  significant  effects  on  both  the  interface  sliding 
stress  and  the  matrix  cracking  stress,  as  well  as  the  ultimate  strength. 

The  purpose  of  this  investigation  is  to  provide  a  3-D  solution  for  the  case  of  a 
matrix  plate  which  is  embedded  with  a  periodic  array  of  fibers  and  is  under  the  action  of  a 
temperature  change.  In  this  study,  a  3-D  analytical  solution  developed  by  Folias  will  be 
specialized  to  the  problem  at  hand  and  will  be  extended  to  account  for  the  presence  of  a 
periodic  extension  of  fibers  by  satisfying  additional  boundary  conditions  within  the  cell 
configuration.  The  solution  is  subsequently  reduced  to  a  matrix  system  with  a  large 
dimension.  Our  previous  experience  has  shown  that  this  type  of  system  is  very  sensitive  to 
even  small  changes  of  the  coefficients.  The  rate  of  convergence,  however,  will  be 
constantly  monitored  by  how  well  the  boundary  conditions  are  satisfied.  Sophisticated 
numerical  algorithms  have  been  developed  to  handle  the  numerical  analysis  which  will  be 
carried  out  in  double  precision.  Thus,  the  3-D  effects  including  those  at  the  boundary  layer 
close  to  the  free  edges  will  be  recovered. 

The  analysis  is  expected  to  provide  us  with  important  information  for  the  temperature 
mismatch.  The  residual  stress  field  along  the  three  different  directions  respectively  (r,  0  and 
z)  will  be  examined  and  the  location  of  maximum  stress  will  be  established  for  the 
subsequent  failure  analysis.  Moreover,  the  relationship  of  the  stresses  versus  different 
ratios  of  the  shear  moduli  (112/M1),  the  fiber  radius  to  separation  distance  between  adjacent 
fibers  (a/b),  the  fiber  radius  to  plate  thickness  (a/h),  and  the  thermal  expansion  coefficients 
(  a(2)/a(i) )  will  be  determined.  Such  information  will  not  only  help  us  understand  the 
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failure  mechanism  better  but  it  will  also  reveal  how  the  critical  stress  to  failure  is  dependant 
upon  the  respective  material  properties. 


CHAPTER  2 


FORMULATION  OF  THE  PROBLEMS 


2.1  One  Fiber  Model 

Consider  the  equilibrium  of  a  body  that  occupies  the  space  I  x  I  <*»,  |  y  |  <»,  |  z  I  ^  and 
contains  two  different  elastic  materials  (fibers  and  matrix)  whose  common  boundary  is 
given  by  a  through-the  thickness  cylindrical  surface  of  radius  r=a  whose  generators  are 
parallel  to  the  z-axis  (see  Figure  2.1).  Both  the  matrix  and  the  fibers  are  considered  to  be 
made  of  homogeneous,  isotropic  and  linearly  elastic  materials.  At  the  interface  perfect 
bonding  is  assumed  to  prevail.  The  surfaces  I  z  I  =h  for  both  regions  are  free  of  stresses 
and  constraints. 

Considering  the  case  of  a  temperature  change,  which  is  the  same  at  any  point  of  the 
body,  and  considering  the  absence  of  body  forces,  the  coupled  differential  equations 
(Navier's  equations)  governing  the  displacement  functions  u(‘) ,  v(') ,  w(>)  ( i=  1, 2,  where 
i=  1  applies  to  the  region  r>a  and  i=  2  to  the  region  r<a)  are: 


■”1  a  ,(') 

m.  -  2  dx 
1 

9  Ji) 

m^  -  2  9y 

3  Ji) 
m.  -  2  9z 


+  V^u^‘^  =  0 

+  =  0 

+  VV‘>:=0 


(2.1) 

(2.2) 

(2.3) 


where 
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2  H  H 

V  s  -2—  +  -3—  +  -2—  (Laplacian  operator) 
3y^  9zr 


Vj  (i=l,2)  =  Poisson's  ratio  of  matrix,  fiber 


Ji)  =  du^‘^ 
®  "  ax 


(i=l,2) 


(2.4) 

(2.5) 

(2.6) 

(2.7) 


The  stress-displacement  relations  including  the  thermal  stresses  are  given  by  Hooke's  law 
as: 


=  2\i.  ( ^  2^i.  AT 

“  ax  m.  -  2  l-2v.  “ 


Cfyy  =  2Hj  ( 


av<^>  .  e^‘>  ,  1+v. 


^ -•  -^)  -  tV- 211:  a  at 

ay  mj-2  l-2v.  “ 


I  1 


ay  ^  ax  ^ 
V  dz  dy  ^ 

”  ax  az  ^ 


(2.8) 

(2.9) 

(2.10) 

(2.11) 

(2.12) 

(2.13) 


where  iii ,  a(‘)  (i=l,2)  are ,  respectively,  the  shear  moduli  and  the  thermal  expansion 
coefficients  of  the  matrix  and  the  fiber. 

As  to  the  boundary  conditions,  one  must  require  that 

asixi  — >oo: 
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o 

II 

Sh 

II 

II 

(2.14) 

aslyl  — »«»; 

(2.15) 

at '  z  I  =  h: 

(2.16) 

atr=a : 

(2.17) 

M)  _ 

^le 

(2.18) 

(2.19) 

(2.20) 

(1)  (2) 

(2.21) 

(1)  (2) 
w  =  w 

(2.22) 

also  the  following  continuity  condition  must  hold: 

at  r=0 :  all  stresses  and  displacements  for  the  inclusion 
must  be  bounded 

(2.23) 

It  will  be  convenient  to  seek  the  solution  in  the  form 

u(0  =  u^KO  +  uW<i) 

(2.24) 

(2.25) 

^(0  =  ^(pXO  ^  ^(c)(i) 

(2.26) 

where  the  component  with  the  superscript  (p)  represents  the  particular  solution  whereas  the 
component  with  the  superscript  (c)  represents  the  complementary  solution,  which  is  a 
function  of  the  x,  y,  and  z  coordinates. 
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The  particular  solution  is  relatively  easy  to  construct.  In  cylindrical  coordinates,  the 
particular  solution  of  the  matrix  (i=l)  is  of  the  form 


"  e- 

(2.27) 

00 

(2.28) 

-(P)(l)  _  -(P)(l)  _  JP)(1)  -  _  A 

10  ra  *’02 

(2.29) 

Integration  of  equations  (2.27)-(2.29)  gives: 

(2.30) 

uf^>  =  0 

(2.31) 

(2.32) 

For  the  fiber  (i=2),  the  particular  solution  has  the  same  form  as  the  matrix.  Thus: 

o«®  =  2b. 

(2.33) 

(2.34) 

-(p)(2)  _  sm  _  -(p)(^  »  „(p)(?>  _  A 

(2.35) 

b„r  +  o®4Tr 
^  \i-2  l+V^  ® 

(2.36) 

(2.37) 
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w 


(P)(2)  _  1 


2v^ 


\^2  1+V. 


■b^z  +  a^^^ATz 


(2.38) 


where  Ao  and  bo  are  constants  to  be  determined  later  from  the  boundary  conditions.  Note 
that  the  above  particular  solution  for  the  fiber  (i=2)  satisfies  the  continuity  conditions  at 
r=0  ,  i.e.,  equation  (2.23). 

In  view  of  the  particular  solution,  one  needs  to  find  six  functions;  v<‘=)(>), 
w<®Xi),(i=i,2)  such  that  they  satisfy  simultaneously  the  partial  differential  equations  (2.1)- 
(2.3)  and  the  following  boundary  conditions: 


at  Izl  =  h : 

-(c)(o.jcXd_^(cxo_n 

^xz  ^yz  ^ 

(2.39) 

atr=a : 

II 

(2.40) 

.(C)(1)  .(c)(3_.(p)(2)  .(PKI) 

(2.41) 

-(c)(1)  .  .(c)(?)  _  _(p)(2)  .(pKD 
'’ra  '’ra  ’’iz  ^Tz 

(2.42) 

s 

1 

cs 

II 

1 

X 

(2.43) 

(CXD  (cX2)  (P)(2)  (pXD 

“e  ■  “0  ~  %  ■  “e 

(2.44) 

w(')(1)-w(c)(2)  =  w(P^2).^(pK1) 

(2.45) 

asr  — >  oo  : 

all  complementary  displacements  and  stresses  fw  the  matrix 
must  vanish.  (2.46) 

atr=0 : 


all  complementary  displacements  and  stresses  for  die  fiber 
must  be  bounded.  (2.47) 
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2.2  Periodic  Fibers  Model 

The  above  discussions  are  based  on  the  assumption  that  the  matrix  plate  is  embedded 
with  one  fiber  only  and  that  the  entire  system  is  subjected  to  a  temperature  change.  Such  a 
model  is  often  used  to  illustrate  how  the  residual  stresses  between  a  fiber  and  a  matrix 
develop  as  a  result  of  thermal  expansion  mismatch.  In  practice,  however,  numerous  fibers 
are  embedded  into  a  matrix  plate.  One  may  conjecture,  therefcnc,  that  the  residual  stress  in  a 
composite  plate  will  depend  strongly  on  the  properties  of  the  fiber/matrix  interface  as  well 
as  on  the  fiber  packing  sequence.  Consequently,  this  model  must  account  for  the  presence 
of  a  periodic  array  of  fibers  that  are  embedded  into  the  matrix  plate.  This  model  is  shown  in 
Figure  2.2. 

Most  of  the  formulas  developed  for  the  one  fiber  model  are  still  applicable  to  the  case  of 
periodic  inclusions  with  the  exception  that  additional  boundary  conditions  need  to  be 
satisfied  within  the  cell  configuration.  Symmetric  conditions  may  be  used  to  simplify  the 
problem  reflecting  the  periodic  fibers  model.  From  Figure  2.2 ,  it  is  clear  that  lines  0-B,  O- 
C,  and  B-C  are  symmetric  boundaries.  Thus,  one  may  consider  the  triangle  OBC  and 
concentrate  on  how  to  deal  with  the  continuity  conditions  on  the  fiber/matrix  interface  and 
on  the  symmetric  boundary  line  B-C. 

For  convenience,  all  displacements  and  stresses  will  be  expressed  in  a  form  that 
automatically  satisfy  the  symmetric  boundary  conditions  along  the  lines  O-B  and  O-C.  It 
remains,  therefore,  for  us  to  satisfy  the  symmetric  boundary  conditions  along  the  line  B-C. 
In  particular. 


c 

n  J  ** 


(2.48) 

(2.49) 


(2.50) 
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Un=  constant  (2  51) 

where  n  and  t  are  the  normal  and  parallel  directions  of  the  line  B-C,  respectively. 

After  separating  the  solution  into  a  complementary  arxl  a  particular  part,  one  may  write; 


1 

II 

(2.52) 

-(cXl)  _  .  -(0(1) 

(2.53) 

C  C 

«(1)  ^  f  ds  =  .  f  ofB  ds  =  -  if” 

n  1  n  1  ^ 

w  ir 

B  B 

(2.54) 

(2.55) 

CHAPTER  3 


METHOD  OF  SOLUTION 


3.1  Conrolementaiv  Solution 

A  general  method  for  constructing  solutions  to  some  3~D  mixed  boundary-value 
problems  that  arise  in  elastostatics  was  developed  by  Folias.  Based  upon  these  results,  one 
may  assume  the  general  form  of  the  solution  to  system  (2.1)-(2.3)  which  automatically 
satisfles  the  boundary  ctxididons  at  the  plate  faces,  equation  (2.39),  to  be  of  the  form: 


Ott 

^  S I  fl®.  +  “i  1 

>  V=1 


X  (-1)"^^  cos(a„  z)  +  y 


(0  dX. 


2  (0 

1  .29X, 


dy 


dx  ntj+l  3x0y 


(3.1) 


(0 


“*1"^  v=i 

.(0 


ndn;"  .  ^  3m- 1  (0  (0 

"  cos(a„z)  +  --^X3  +X2 


n=l 


ax 


(® 


dx 

2  (0 

.  ^3 


dy  +1  dx 


••  .. 
™>  ^  v=l 


(3.2) 
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m.+l  By 


Furthermore,  the  stresses  are  given  by  equations  (2.8)-(2.13).  These  expressions, 
however,  do  not  include  the  thennal  stress  con^nent* 


=  nf  f,(P, I) +^42(111-1) f,(p. 2) +  m.fj(P.z)]  ) 


m.-2  ^ 

»  V=1 


an® 


(9  2  (9 

aXj  a  X3 


“  dxdy  dx  dx 

,  2  <,  1 

nij+l  ay  nij+l  ax^ay 


s:  ®w  °  X  *  ^1  >^'“1'*'  'i<P'  ^>1  > 


m.-2 

1  v=l 


(9  (9 


V/ixn3H„  ,  ,  2m  aXj  aXj 

^  ayZ  m+1  ay  ax 


2  (9  3  (9 

ax^  mj+1  ax^dy 


^_(C)(9_  frft 


»  aTr'  ^ 

Ijr  “  m?T  I^ay  ) 

•  i  v=l  ^ 

^ ,  ,.n ,  a^l'*  I  2  u®  1  ,  “r’ 
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(i)  2  (i)  3  (») 

,9^2  1 

3x  ^  dxdy  nij+l 


=  (  f,(P.z).f/P.z) ) 

t  V=1  ” 

+  i£(-l)"^sm(a„« 


1  v=l 

.„(«■) 

-  y  X 

Ttel  ^ 


where 


an  =  x  ’  "=^’2’^ . . 


(3.10) 


andPv  are  the  roots  of  the  equation 


sin(2pvh)  =  -(2p^h) 


(3.11) 


Hv<>)  and  Hn(>)  are  functions  of  x  and  y  (  or  r  and  6  depending  upon  the  system  of 
coordinates  used)  that  satisfy  the  reduced  wave  equation: 


7)^  2  (0 

-Pv)Hv  =0 

dx*^  dy^ 


(3.12) 


,2  2 


(:r2  +  #T-“5)Hl=0 


9x^  3y 


(3.13) 
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Xi(‘) ,  X2(*) ,  and  XsCO  arc  functions  of  x  and  y  ( or  r  and  0  depending  upon  the  system  of 
coordinates  used)  that  satisfy  Laplace's  equation: 


dx  3y 


In  addition,  Xi(»)  and  X2<‘)  satisfy  the  relation: 


(3.14) 


dx  By 

Finally, 

f  j(Pv  z)  =  cos(Pv  h)  cos(Pv  z) 

f^CPv  z)  3  Pv  h  sin(Py  h)  cos(Pv  z)  -  Pv  z  cos(pv  h)  sin(Pv  z) 
fjCPv  z)  3  cos(pv  h)  sin(pv  z) 

f^  (Pv  z)  3  py  h  sin(pv  h)  sin(Py  z)  +  Py  z  cos(Pv  h)  cos(pv  z) 


(3.15) 

(3.16) 

(3.17) 

(3.18) 

(3.19) 


Because  the  above  complementary  solution  satisfies  the  boundary  conditions  (2.39),  all 
that  remains  now  is  to  satisfy  the  continuity  conditions  at  the  tiber/matrix  interface, 
equations  (2.40)-(2.45),  and  for  the  periodic  fibers  model  only,  the  additional  symmetric 
boundary  conditions  along  the  line  B-C,  equations  (2.S2)-(2.5S),  as  well  as  continuity 
conditions  (2.46)-(2.47). 


3.2  One  Fiber  Model 


In  the  case  of  one  fiber  model,  there  is  only  one  kind  of  boundary  condition  that  needs 
to  be  satisfied,  that  is  the  continuity  condition  at  r=a  .  By  using  the  coordinate 
transformations  of  displacements  and  stresses,  the  complementary  solution  can  be 
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transformed  firom  the  rectangular  to  cylindrical  coordinates  in  order  to  take  advantage  of  the 
symmetry.  Thus,  equations  (2.40)-(2.45)  may  be  written  in  the  following  form : 

sin%  )  +  cos^  )  +  sin(2e)(T<;>(« 

Jc)(2)  s 
■>  ^ 

A„ 

=2b,-^ 

(3.20) 

1  sin(20)  ^  sin(28)  +  cos(2e)  -  x^^^ 

=0 

(3.21) 

sine  (X<«‘>  -  )  +  cose  (T«»  -  »  )  =  0 

(3.22) 

sine  (u^X”  -  n«® )  +  cose  (v«<'>  -  ) 

(3.23) 

cose  (u«<'>  -  „«® )  -  sine  ( v“>'»  -  v‘'xa  > = o 

(3.24) 

^,0X.) .  „(CK2)  ^1  2v,  ^  ^ 

M-2  I+V2  ” 

(3.25) 

Taking  into  account  the  symmetry  of  the  problem  ( 0  independence ),  one  finds: 

h"’  =  “o. 

(3.26) 

Hf  =  b^yM 

(3.27) 

(3.28) 
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=  (3.29) 

and  all  X=0,  where  Iq  and  Ko  are,  respectively,  the  modified  Bessel  functions  of  the  first 
and  second  kind  of  order  zero;  aov  and  bov  are  complex  constants;  and  con ,  don  are  real 
constants. 

It  is  found  convenient  to  nondimensionalize  the  expressions  by  introducing  the 
following  definitions: 


(3.30) 

^Ov  ■  ^(1)  ^Ov  Pv  lo^Pva) 

(3.31) 

Bon=  (1? 

AT 

(3.32) 

®0n  ”  (1)  ‘^On  ^ 

OC  M 

(3.33) 

A*  -  1  A 

0  ~  (1)  9  0 
^  AT  Pj  a^ 

(3.34) 

^0  ~  (1)  ^0 
^  AT  ^2 

(3.35) 

Substituting  equations  (3.26)-(3.35)  into  (3.20)-(3.25),  one  has 


1  (1)  l^n(3v^) 

zL  aL  I  2f,(M  +  I2(ni,  ■  l)f,(M  +  ni,f2<P.z)l  ^.2 - ) 

1  (2)  1^'? 

■  SJ-  X  V  1  2f,(P.^)  +  +  n^yM)  I  ^ 

2  v=l  Hv  fofPv^) 
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^T-tTr^o-® 


(3.36) 


®On  t  J-  2  ^  ^  cos(ce„z) 


2  ot^a 


(3.37) 


m,  ^  (1)  K»(Pv^) 

^  S  ^ov  KgOva)  ^  ^ 


T,  I  f3(PvZ)  +  f4(PvZ)  ]  ^  =  0 

^-2  ^  Pv  lo(Pva)  ^  ^ 


(3.38) 


1  X*'  (1)  1  ^n(Pva) 

^  X  A,  [  2(m,.  I)f,(M  .  ™,f,(M  1 

1  (2)  1  TnCPya) 

A,  (2<^-  „f.(M^«.,f,(M]^^ 


4.^,;.,.4=o 

2  1+v  0  (1) 

2  a 


p;a  V 


(3.39) 


litl  '"’i"  -  “ifjfM  1 

ir=l 


(3.40) 
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1 

1112-2 


S  ^^“2’  -  ni2f4(PvZ)  ] 

v^l  Pv® 


a 


(2) 


a 


(1) 


3.3  Periodic  Fibers  Model 


(3.41) 


For  the  periodic  fibers  model,  one  needs  to  consider  the  additional  boundary  condition 
along  the  line  B-C.  From  Figure  2.2,  it  is  noted  that  the  displacement  and  stress 
distributions  are  now  0-dependent.  In  order  to  satisfy  the  symmetric  conditions  at  6=0  and 
0=7c/2  respectively  (see  Fig.2.2),  one  needs  to  chose  the  solution  in  the  form  of  a 
summation  of  cos(2k6)  or  sin(2k6).  Thus,  the  symmetric  conditions  along  lines  O-B  and 
O-C  will  be  satisfied  automatically.  Moreover,  the  summation  provides  a  sufficient  number 
of  coefficients  so  that  the  symmetric  condition  along  the  line  B-C  may  be  satisfied.  Thus, 
the  complementary  solution  may  be  chosen  in  the  form: 


h1«  = 

X  Sy  cos(2ke) 

ksO 

(3.42) 

=  Z  ^kv  y  M  cos(2ke) 

(3.43) 

h1«= 

^  ^kn  )  sin(2k0) 

fc=0 

(3.44) 

sin(2k0) 

(3.45) 

^  ^  Sin[(2k+l)01  +  ^  sin[(2k-l)0] 

k=l  ^  Ic=l 

(3.46) 
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C  =  X  J?  1  cos[(2k+l)e]  -  2^  cos[(2k  -1)0] 

(3.47) 

4‘’  =  ^:^«>st(2k+l)01 

(3.48) 

^  sin[(2k-l)0] 

k=l 

(3.49) 

=  -  ^h^T2^-^cos[(2k-l)e] 
k=l 

(3.50) 

=  ^ikr^^^cos[(2k+l)0] 
k=0 

(3.51) 

where  Iik  and  Kik  are ,  respectively,  the  modified  Bessel  functions  of  the  first  and  second 
kind  of  order  2k,  and  akv,  bkv,  Ckn>  dkn  (  k=0,l,2...M;  v,n=l,2,3....).  Ck  ,  fk  ,  gk  .  hk  ,  ik 
are  arfeitrary  constants  that  are  to  be  detennined  later. 

For  numerical  calculations,  it  is  found  convenient  at  this  stage  to  nondimensionalize  the 
expressions  by  introducing  the  following  definitions: 


(1)  _  ^  P'' 

a  AT 

(3.52) 

a  AT 

(3.53) 

(3.54) 

®ln= 

a  AT 

(3.55) 

(3.56) 
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>  fl2k-2 

f*  _  a  r 

(3.57) 

(3.58) 

*  a2k-2 

^k  ~  (1)  \ 

^  AT 

(3.59) 

.*  a^*^  • 

“  (1)  *k 

a  AT 

(3.60) 

Now  the  solutions  of  both  particular  and  complementary  parts  have  been  found.  The 
only  problem  remaining  is  to  determine  the  unknown  coefficients.  Thus  one  needs  to 
concentrate  on  how  to  deal  with  both  the  continuity  conditions  at  the  fiber/matrix  interface 
(r=a)  and  the  symmetric  boundary  conditions  along  the  line  B-C  (see  Fig.2.2)  For  the  first 
boundary,  one  can  use  a  Fourier  expansion  in  order  to  eliminate  the  z-dependency  and  thus 
recover  a  set  of  simultaneous  linear  equations  of  the  unknown  coefficients.  For  the  second 
boundary  (line  B-C),  one  can  use  an  appropriate  numerical  scheme,  which  may  match  the 
symmetric  condition. 


CHAPTER  4 


NUMERICAL  CALCULATIONS 


4.1  Fourier  Expansions  for  One  Fiber  Model 
In  order  to  simplify  the  numerical  calculations,  (xie  must  first  eliminate  the  z-dependence 
at  the  boundary  r=a  and  thus  recover  a  set  of  simultaneous  linear  equations  with  the 
coefficients  as  the  unknowns.  It  is  quite  interesting  to  note  that  Bon(^)  and  Bon^^)  are  not 
coupled  in  the  system  (3.36)-(3.41).  That  means  equations  (3.37)  and  (3.40)  may  be 
treated  separately  from  the  remaining  four  boundary  conditions  and  the  following  solution 
is  obvious: 


^(1)  ^(2)  ^ 


(4.1) 


In  view  of  equation  (4.1),  two  of  the  boundary  conditions  are  now  automatically 
satisfied  and  only  four  are  remaining.  The  relative  balance  between  unknowns  and  number 
of  equations  does  not  change  because  two  group  of  constants  are  satisfied  by  (4.1)  and 
these  two  relative  group  of  equatitxis  are  autcmiatically  satisfied. 

The  technique  for  solving  such  a  system  is  given  by  Folias.  More  specifically,  one 
expands  the  functions  fi(^  z),  f2(^  z),  and  (r/a)^  in  terms  of  cos(an  z)  and  the  function 
f3(Pv  z),  f4(^  z)  and  (z/a)  in  terms  of  sinfOn  z).  Equating  next  the  coefficients  of  similar 
functions  of  z ,  one  arrives  at  the  following  equations,  which  may  then  be  used  directly  for 
the  numerical  calculations.  The  reader  may  note  that  for  convenience,  the  following 
notati(»i  is  introduced : 
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r  =  ■  ^  • 

Py-«5 


™1  *^  1 
1  V=1 


^0  .  u*  ^^2 

T-^^oTf; 


2  'V  .®  ■  ,  .  'o<P»*)  .  1*2 

inr 


uXcM)  P^(M)  ’‘■ 


2  ^  2  ^  1  yPv^^ 

Ko<P.» '  "“p^VP.*) 


•ft 

Ia^ 


_  '^O  o£^..  J 

“T"  1+v  0  0) 

Z  l+Vj  QjV 


1  v=l  Pv^ 


*  v=l  Pv^ 


1 

=  -i— (l..Si — ) 
Ona^  „(!)' 


,  ^  m  k;;(M 

Si  A^  r„,  {-l-K  -in,+l+m,r„  )  -y5 - ) 

>1-2  PX(P,a) 


(4.6) 
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1 

m2-2 


v=l  Pv 


fpCPya) 

Io(Pva) 


m. 


mj-2 


V  “n  Ko(p^a) 
^  0vp2  "VK„(p^a) 


_J_Y  A^^^fbLT^  yPv^>^^2_o 

"^-2^  "''lo(Pva)h" 


(4.7) 


(4.8) 


This  linear  system  can  be  solved  by  a  modified  Gaussian  elimination  method  with 
pivoting.  Double  precision  is  used  throughout  the  computational  analysis  because  of  the 
high  sensitivity  of  the  system. 

The  reader  can  show  that  the  solution,  by  virtue  of  its  construction,  satisfies  (i)  Navier’s 
equations  (2.1)'(2.3),  (ii)  all  the  boundary  conditions  and  (iii)  all  the  continuity  conditions. 
Moreover,  the  following  three  special  cases  will  serve  as  a  limit  check  of  the  numerical 
results.  The  three  special  limit  cases  are  a  continuous  plate;  a  thin  plate;  and  a  plate  with  a 
cylindrical  hole. 

(1)  Continuous  plate.  If  vi  =  V2 ,  Pi  =  P2  •  Then: 


(4.9) 


So  the  solutions  reduce  to  the  particular  solutions  and  all  of  the  residual  stresses  go  to  zero 
and  the  displacements  reduce  to  the  linear  forms  along  the  different  directions. 

(2)  Thin  plate.  If  h/a  «  0,  then  the  right  hand  side  of  equation  (4.6)  will  be; 


iL  =  XiL_»o 

a„a  ~  aji  a  nn  a 


as(^)-^0 


(4.10) 


from  equations  (4.3)-(4.8),  one  has  : 
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.(1)  .(2)  (1)  (2) 
A(v=Ao.=®(h=®(h=0 


(4.11) 


It  is  seen  that  (me  retmvers  precisely  the  particular  solution,  that  is  the  plane  stress  solution 
given  by  Delale. 

(31  Ot^lindrical  hole.  In  the  case  of  which  is  numerically  represented  by  a  very 

very  soft  inclusion,  i.e.,  1=0.00001.  One  can  see  from  equations  (4.7)-(4.8)  and  from 

equation  (4.3)  that: 


a;=o 


(4.12) 

(4.13) 


which  suggests,  therefore,  that  in  the  case  of  a  matrix  plate  with  a  cylindrical  hole  there  is 
no  residual  stress  present  as  a  result  of  a  change  in  temperature. 


As  for  pericxlic  fibers  nKxlel,  there  are  two  di^erent  kinds  of  boundary  conditions  that 
need  to  be  satisfied.  Similarly,  for  the  continuity  conditions  at  r=a,  Fourier  expansions  are 
used  to  recover  sets  of  linear  equations  with  unknown  coefficients.  Without  going  into 
mathematical  details,  the  following  equations  can  be  used  as  the  boundary  conditions  at 


2  •y  1  _2_^ 

VM  ■  ">2-2  ^  ‘^B?.VP.» 


2^  mj+1  2^m2+l  ^  « 


.^l/3m2-l 


2  l+Vj  «  „») 


(4.14) 
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P.Ko(3.a)  ^'^*=1  P.Io<P.a) 


r^tM)  ,m 


m,-2  ^ 

1  V=1 


(°^i'^  -hDe*  +i*i^ 

^+1  ^®0  o^i  2  onj 


(4.15) 


1  ^  (1)  r  K' Oya) 

X  V  rT  ^  -m  +l+m  r„y )  - 
^  Pva  KQ(Pva)Py 


n-2^"Ov  p  '  "1  “1 

1  V=1 

2  V=1  PyS 


r„(P.a) 

lo<P.a) 


(4.16) 


ni^-2  X  V  2"  (™r^+™l^nv)rn 
1  V=1  Pya 


■"  i  ^  ^  (mj-i+mjr„)  r „ 

^  V=1  Pya 


_2 _ L_i*+_J_(i.^) 

Onara^+l  0  o^a^  ^(i)  ^ 


(4.17) 


1  (1)  Kn(Pv^) 

I'nv  (  -1+  ( -mj+l+m^r^y)  ) 

pX(M) 


1  ’'f'  (2)  r  ^n(Pva) 

;^XAo.^«-l4(-n,.Hn,r.)^)^  =  0 


(4.18) 


V  ^  KQ(Pya) 
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1  V  A®  “nyJ 


(4.19) 


2  'V  1  •^k<Pv4>  2  V  a!^’  ’ 

”72  ^  ^  pja  ■  ”2-2  ^  pja  ‘aCM) 


1  ,2m,-* 


I  ,2”,-* 


i*^^2L.l)gI.i(^.2k.l)gJ.,.^i)^^2k..K2k,g;,, 


*  .  l/h^2  1 


(4.20) 


i_  Y  A5'^__k_+_4_  Y 


,<2)  k 

kv  .k2 


= ■  i  i<- iff  ■  T 


I  3m, -1 


1  /hN2  1 


^  i<$r  ^  i<-  ^  -  ^  -’>  *^1  ■  ^ 


.•.1,  2”2-1 


•  l/h^2  1 


(4.21) 


1  Y  f  2  ^2k(M  >  _J_Ya^^^2  2  ^2k(M  .^2 

^  ^  K2.(M  ^  ■  in,-2  ^  'V  t-  -  p2  I^(p^a)  V, 


=  (2k+l)(^J^+k+l)  g;+(2k-l)(k+l)  g;  ,  +  -|<^)^-j^2k-l)(2k)(2k+l)  g;., 


•  1^2 


+  (2k+l)(l-k)  i;^+  (2k-l)(^-k+l)  i;,^ 


+ 1  (^)^^j^2k+l)(2k-l)(2k)  -  (l-2k) 


+  (l+2k)cj+(l-2k)h;;^ 


(4.22) 


4  Y  ^1)  ,  -k  ,  k  I 

pJa^2k(Pva)  (p^a)^ 


(2)  f  -k  ^2k(Pv^)  ,  k  ,  ^^2 

^  l  V  '  <>  '  II 


'2  “  v=l 


l2k^Pva)  (p^a) 


,2  ' 


=  (2k+i)(^  +k+i)  g*  +(2k-i)k  g;; , + (2k-i)(2k)(2k+i)  g*. 

m,+l  K-i  3  a  m.+l  * 


+(2k+i)ki;^  -(2k-i)(^-w)i;.,^ 


-  T<T)';;r^T  (2k+l)(2k-l)(2k)  -  (2k-l)  fj 


3  'a'  m^+l 


+  (2k+l)c;;+(2k-l)h;^ 


(4.23) 


--—  ^  A.  — ^  {  -in,+l+in,r„y  ) 
m,-2  kv  2  ‘  1  1  ’ 

*  V=1  Pya 


K2^(Pya) 

K2^(M) 


2_'f'  a^L"  (  „  *14.™  r  1 

1.-2  2-  V  „2  >  I 


m,-2  ^ 

2  v=l 


B;„^^^^(2k-l)(2k)  g;, 

(a^)  (a^)  (a„a) 


2  1 
(a„a)^ 


(2k+l)(2k)C  =  0 


(4.24) 
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-8 


m,-2 


(1)  k 


1  v=i  (M) 


,0^.2  t‘™2‘*'^'*’®2^nv}rnv 


{  -m^+I+mjrnv  )r„. 


“2-^  ^  (M) 


(i)_i_j^^M  pO  1 

'"  KaCV)  ■  *"  aja  W) 

+  -l^-L(2k-l)(2k)  g^,  +  _l_^(2k-l)(2k)  ■;= 


0^1-2  —  -  o‘„ 

1  V=1  Pv^ 


PCa 

1  1 


m,-2  — 
1  V=1 


^  (1)  J^2u(Pv^) 

y  aV  (-l+[-m,+l+m,r^  1  - ) 

^  P:  Ka(P,a) 


Y  A®  %  ( .IH-nu.Um,r„  ]  ,  h 

(1),  k  k  ,  (^  k  ,y“na)  1  ^l; 

■"  aja  K2,(a„a)  ^2^  ti^(a„a)  ‘  a  V 


2  1 


l^(2k-l)(2k)(2k+l)g;, 

(a„a)"™i+l 


^ (2k+l)(2k-l)(2k)  i!  ^^2  _  0 


(4.25) 


(4.26) 


(4.27) 
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rm+1+mr  1  .-21c  .  2k  , 

I  9  ^  '  ^kv  ^nv  I  roj+l+nijr„y  ]  I  /Q  ^  2  ^ 

pJa*^2k(Pva)  (p^a)^ 


m,-2  — 

1  v=l 


A^’  r„.  [-m3+l+  ]  (^- 

•  v=i  Bya  ^ 


r2k^2k(Pva)^  2k 


P^a  ^2k^Pva)  (Pya)^ 


(1)  _j_K^^  4k^ 

■  ■"  aja  K^(a„a)  ■  (a„a)^  ■  2  ' 

,  p(^  ,  1  ^'ac^“n^)  4k^  .1}!^ 

aja  l2k(«na)  (a^a)^  2 

(a„a) 


+-^-;;;^(2k+i)(2k-i)(2k)i!^  =o 

(a„a)^™2+^ 


(4.28) 


illlLV  ^2k^Pv^)  .  V  «n 

"%k(M 


(pya)  )l2 
.(Pya) 


+  B 


k  p<^  k  ^2  _Q 

kn  a„a  ■  to  a„a  Hj 


(4.29) 


As  for  another  group  of  boundary  conditions  at  line  B-C,  Fourier  expansions  are  also  used 
to  eliiiunate  the  z-dependency.  Then  the  collocation  method  is  used  to  force  the  conditions 
(2.52)-(2.55)  to  be  satisfied.  Without  going  into  mathematical  details,  the  following 
equations  can  be  used  directly  in  the  numerical  calculations  where  the  following  definitions 
have  been  used: 


b 

sin(0p+cos(0j) 


(4.30) 
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n 


^  ^  number  of  points  along  B-C  line) 


5S  -  _ sin(6e) 

‘  ^  sin(2e.+8e)+cos(8e) 


(4.31) 


(4.32) 


(4.33) 


at  B-C  line: 


Mo, 

Kjj(P.a)pJ  pj^.  Ka(M 

4k(k+l)  K  (P,4^ 


M  „ 

^2k^Pv^i) 

'  Kj^(P.a)  pj  pj|.  K2k(P.a) 

4k(k-l)  J^ac^Py^j)  , 

’o.5,)'  K3.(P^)  l««"»-2)e^ 

+  ^  B^^^  f  k-K).5  ^  k(l+2k)  i  ,  K^(a„^;) 

a%  K^(a„a) 

+  'f'  B^*^  f  jt-0.5  k(l-2k)  ]  -  K3^(a„  U 

S  K^aKa)  '  4^  '  “«2k-2)ej 

(2''-l)(2k)(2k+l)(|.)»"2g^jOos((2k+2)9j|  =0 


(4.34) 
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M 


il.w 


^kv  2  ^  ^nv  ( 


m,-2  ^2 

1  k=o  v=i  Pva 


{  cos[(2k+l)0jl  +  sin[(2k+l)0jl  ) 


K^kCPv^  2k*^2k(Pv4i) 
K2^(P,a)  ■  4. 


] 


K^O.a)  4,  K^(M) 


{  cos[(2k-l)0jl  -  sin[(2k-l)0jl  } 


M 


X  X  4“  ^  -mi+l+™irnv)  ^  ^  cosfe 

1  k=o  v=i  Pva  ^acCPv^)  ^ 


+  :S'^b“’-!-(  2k  Ka.(a„4i)  Ka(a„4i) 

“  a  '"  aja  5i 

{  cos[(2k+l)0jl  +  sin[(2k+l)0j]  ) 

b(>)  1  r2k  ^2kK^i)  , 

^1^  «5a  K2,(a„a)  K^^Ca^a) 

(  cos[(2k-l)0^  .  sin[(2k-l)0p  ) 

+  '^  B*”  2k  '^ac<“.^*)  ^,(k7L) 


2J2  1 

(a„a)^  “ 


i 


(2k-l)(2k)(f)^^‘g* 


{cos[(2k+l)0jl  +  sin[(2k+l)0jl  ) 


1 


"^1+1  (a„a)"  - 


^  (2k-l)(2k)(4)^^^g;  jCosfe 


=  0 


(4.35) 
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^  f  t  A*"  r„  (  f  {  . 

pJk^CM  (P,4/  '^3k(M  ‘  ^ 

r..(-n...H„.r^)  it  i 

^  ^  ti  p.K^(M)  PSi 


m,-2  — -  — 

1  k=0  v=l 


7  '"‘*'^*2  ^  1  »Si  sta((2k+2)e^ ) 

^2k^Pv^)  (P^  K2|.(Pva) 


M 


^  X  £  ^kv  r„v(  -m,+l+m,r„,)  {  £  [  - 
1  k=0  v=l 


K"^(Pv^i) 


"  p,  K^^CPva) 


,  l-4k  ^ac^Py 


4k(k-l)  K2k(fiv^i) 


]  SS.  sin[(2k-2)0.J ) 


fc=0  i=l 


r  k.K).5  '^'ac<“n^i>  t(U2t)  Ka(«„^ 


+ 


+ 


4  K^(a„a) 


]  6Sj  sin[(2k+2)ejl  ) 


4  K^(a„a) 


]  SS.  sin[(2k-2)ejl  } 


1 


(a„a)^  fc=i 


1 


(2k-l)(2kX2k+l)  gj 


i=l 


V 
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sin[(2k+2)ejl  }  =0 


(4.36) 


2 -1-2  6  "  pjK^(Pva)  ^2k(3va) 


4k(k+l)  Ka^(Pv  ^j) 
^2k(Pv4) 


)  cos[(2k+2)ejl 


1  k=0  v=l 


(1)  -  ^  |.4]^  K^j^Cpv^j) 

p5lCjj(p,a) 


4k(t-l)  l^ac(P.^i) 

(P,4/  '^a<M 


}  cos[(2k-2)6.]  ) 


(2k+l)(^l^+k+l)  cos[(2k+2)0jl 

1 


+  ^  (2k-l)k  if-f  g'  j  cos[(2k+2)ejl 

1^1  Sj 

+ S  (2k-l)(2k)(2k+l)  cos[(2k+2)0jl 

1  1^1 

+  (2k+l)(^)^^S;;  cos[(2k+2)0jl 


+  ^  (2k-l)  cos[(2k-2)0^ 


iA;(|-)^  cos(20p  =  O 


(4.37) 
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72  1  y  Y  ,»>  1  .  2k  *^2k<Pv^i>  ■ 

2  '"pja  Kj^CM) 

{  cos[(2k+l)ej  +  sint(2k+l)0j| ) 

,72  1  f  t  .-l'”  »  I  i  ^  I 

2  “r^  ^  pja  5,  ^(.M 

(  cost(2k-l)e,l  -  sui[(2k-l)e,| ) 


-2-'^  y  A"’-!-^^^^cos(te 


+  ^  y  1  +2k+l]  (A)*^'  g;  (cos[(2k+l)6J  +  sii.l(2k+l)e,l  ) 


+  ^  y  (2k-l)  (f)’‘‘‘g;.,  (cos[(2k+l)0|)  +  sm[(2k+l)0,|  ) 

^  k=l 

*  '  ^-1  >  «>s[<2k-l)0il  -  sin[(2k-l)0J  ) 

15=1  1 

■"  ^  (2k-l)(2k)  g;,{cos[(2k+l)0jl+sin[(2k+l)0jl) 

1  fc=l 

fc=0  1  s 


k=l  1  S 


a  x2k-l 


J  (2k-lX2k)(i)»^'  g;,,cos(f  ) 
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-  ^  (2k-l)  ^  (^)^5SjCos(2ke^  } 

k!=l  tl 

-  ^  (2k+l)(  ^  +k+l)  gj  1  ^  8S,  sin[(2k+2)ej  ) 

fc=0  1  '•i 

-^(2k-l)k^.j  sin[(2k+2)ejl  } 

k!=l  i=l 

(2k-l)(2k)(2k+l)g;,  (  ^  (|-)*^^6SjSin[(2k+2)ejl  } 

1  k=l  i=l 

-  ^  (2k+l)  ej  {  ^  sin[(2k+2)ejJ  ) 

fc=l  i=l  S 

+  ^  (2k-l)  ^  4)^'^5Sj  sin[(2k-2)ejl  } 

k=l  i=l 


+  sin(2ej)  =0 

(4.39) 

These  two  linear  systems  (2.14)-(2.39)  can  be  dealt  with  by  a  modified  Gaussian 
elimination  method  with  pivoting.  Also  double  precision  is  used  throughout  the  numerical 
analysis  because  of  the  high  sensitivity  of  the  system. 

4.3  Convergence  of  the  Fourier  Expansions 
The  convergence  of  the  Fourier  expansions  obtained  above  determines  whether  the 
solution  approach  is  a  success  or  a  failure.  In  general,  including  the  problem  under 
consideration,  the  convergence  of  the  solution  is  constantly  noonitored  by  (i)  how  well  the 
series  coefficients  converge  and  (ii)  how  well  the  boundary  conditions  are  satisfied  within 
the  specified  allowable  error. 
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Many  calculations  for  various  parameter  combinations  have  been  completed  in  this 
thesis.  Theoretically  speaking,  the  large  linear  system  (4.14)-(4.39)  can  be  solved  by  the 
use  of  a  general  numerical  scheme  with  pivoting.  However,  because  in  some  systems  the 
rate  of  convergence  may  be  slow,  the  weighted  residual  method  is  also  used  to  eliminate 
possible  ill  conditions. 

A  typical  example  of  the  convergent  coefficients  for  three  different  sets  of  roots  is  given 
in  Table  4.1.  From  the  table,  it  is  observed  that  the  coefficients  converge  slowly  and  that 
200  or  300  roots  are  sufficient  to  provide  us  a  solution  with  the  desired  accuracy.  Figures 
4. 1  and  4.2  show  how  well  the  boundary  conditions  for  the  displacement  u^  and  the  stress 
o„  are  satisfied,  at  the  interface  r=a,  for  three  different  sets  of  roots.  It  is  noted  that  as  the 
number  of  roots  increases,  the  error  on  the  boundary  conditions  decreases  accordingly. 
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Table  4.1  The  convergence  of  coefficients  with  different  roots 


(vi^.34,  v2=0.22,  yi-2/\ii\=\6.61  and  a/h=0.05) 
The  real  pan  of  the  coefficirats 


V 

200  roots  case 

100  roots  case 

50  roots  case 

I 

-1.670072099D-3 

-1.668732365D-3 

-1.664213512D-3 

3 

-6.191873001D-4 

-6.188071283D-4 

-6.163917258D-4 

5 

-3.2473671 13D-4 

-3.24585 1946D-4 

-3.225395566D-4 

7 

-1.963034675D-4 

-l.%2027912D-4 

-1.941973547D-4 

9 

-1.298964396D-4 

-1.297835499D-4 

-1.2771 12265D-4 

11 

-9.167471579D-5 

-9.152634067D-5 

-8.934720024D-5 

13 

-6.7900337 12D-5 

-6.77091 1603D-5 

-6.540784728D-5 

15 

-5.220702888D-5 

-5.197208492D-5 

-4.9542771 14D-5 

17 

-4.1 349622 15D-5 

-4.107248438D-5 

-3.851230089D-5 

Figure  4. 1 .  Error  of  bound 
where  vi=0,34,  V2=0.22, 


Figure  4.2.  Error  of  boundary  condition  at  r=a  for  On 
where  vi=0.34,  V2=0.22,  ^i2/Hi=16.67  and  a/h=0.05 


CHAPTER  5 


RESULTS  AND  DISCUSSION 

5.1  One  Fiber  Model 

Once  the  coefficients  have  numerically  been  determined  firam  equations  (4.3)-(4.8),  the 
displacements  and  stresses  may  then  be  computed  at  any  point  in  the  plate.  In  the  case  of 
one  fiber,  the  displacement  and  stress  fields  arc  independent  of  6  because  of  symmetry. 
Thus,  the  displacement  ue  and  stresses  ,  Xqz  will  be  automatically  zero.  On  the  other 
hand,  because  the  ratio  of  thermal  expansion  coefficients  a(2)/a(i)  presents  a  linear  action 
to  the  displacement  and  stress  field,  the  ratio  a(2)/a(i)=0.075  is  fixed  throughout  the 
discussion  in  this  thesis.  It  is  assumed  in  the  thesis  that  the  residual  stresses  arose  from 
thermal  contraction  mismatch  upon  cooling. 

Because  the  residual  stresses  are  induced  by  the  mismatch  of  the  different  material 
properties  of  the  fiber  and  matrix,  such  as  the  differences  of  the  thermal  expansion 
coefficients,  the  shear  moduli  and  Poisson’s  ratios,  the  interface  of  fiber  and  matrix  is  of 
greatest  practical  interest.  Plots  at  r=a  are  given  from  Fig.5.1  to  Fig.5.18  for  the 
displacements  and  stresses  as  functions  of  z.  In  order  to  examine  how  the  material 
properties  affect  the  field  of  displacement  and  stress,  two  different  ratios  of  shear  moduli 
and  Poisson's  ratios  are  chosen  to  represent  two  practical  composite  materials.  One  is  the 
glass  fiber/epoxy  matrix  composite  with  a  shear  moduli  ratio  of  p.2/lti= 16.67  and 
Poisson's  ratios  of  vi=0.34,  V2=0.22.  The  other  is  boron  fiber/aluminum  matrix  composite 
with  the  shear  moduli  ratio  of  p,2/^i=6.3  and  Poisson's  ratios  vi=0.2,  V2=0.33.  Each  plot 
here  consists  of  three  curves  corresponding  to  different  geometrical  ratios  a/h.  Superscripts 
areusedtodistinguish 
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Figure  5.1.  Displacement  Ur  at  p=a  vs  zDci  for  one  fiber  model 
where  vi=0.34,V2=0.22,p2/|Ai=l6.67 
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Figure  5.2.  Displacement  Ur  at  r=a  vs  z/h  far  one  fiber  model 
where  vi=0.2,V2=0.33,H2/4i=6.3 


Figure  5.3.  Displacement  w  at  r=a  vs  z/h  for  one  fiber  model 


where  vi=0.34,V2=0.22,|i2/lii=16.67 


Figure  5.4.  Displacement  w  at  r=a  vs  z/h  for  one  fiber  model 
where  vi=0.2,V2=0.33,p,2/M-i=6.3 


Figure  5.5.  Stress  On  at  r=a  vs  z/h  for  one  fiber  model 
where  vj  =0. 34,V2=0. 22,^2^41 = 1 6, 67 


Figure  5.7.  Stress  aeeCi)  at  r=a  vs  ?/h  for  one  fiber  model 


where  vi=0.34,V2=0.22,^2/M-i=16.67 


Figure  5.8.  Stress  aee(i)  at  r=a  vs  zlYi  for  one  fiber  model 


where  vi=0.2,V2=0.33,Li2/ui=6.3 


Figure  5.9.  Stress  azz(i)  at  r=a  vs  z/h  for  one  fiber  model 
where  vi=0.34,V2=0.22,ii2/M-i=16.67 


Figure  5.10.  Stress  OzzO)  at  r=a  vs  z/h  for  one  fiber  model 


where  vi=0.2,V2=0.33,U2/Ui=6.3 


Figure  5.11.  Stress  Xn  at  r=a  vs  z/h  for  one  fiber  model 
where  vj  =0.34, V2=0.22,^2/U]=1  6.67 


Figure  5.13.  Stress  Ozz^^)  at  r=a  vs  z/h  for  one  fiber  model 
where  vi=0.34,V2=0.22,^L2/M^i= 16.67 


Figure  5.14.  Stress  at  r=a  vs  z/h  for  one  fiber  model 
where  vi=0.2,V2=0.33,H2/lti=6.3 


Figure  5. 16.  The  ratio  of  stress  (3D/2D)  at  r=a  vs  z/h  for  one  fiber  model 

where  vi=0.2,V2=0.33,H2/|ii=6.3 


Figure  5.17.  The  ratio  of  stress  CoeOy  (3D/2D)  at  i^a  vs  z/h  for  one  fiber  model 

where  vi  =0. 34, V2=0. 22,p.2/Hi=16.67 
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whether  the  displacement  or  stress  is  in  the  fiber  or  in  the  matrix.  If  the  displacement  or  the 
stress  is  the  same  at  the  interface  r=a,  then  the  superscripts  will  be  omitted. 

Figs.  5. 1-5.2  show  the  displacement  Ur  vs  z/h  for  two  different  composites.  The 
positive  value  of  u,  in  most  of  the  center  region  indicates  that  as  the  plate  cools  down,  both 
fiber  and  matrix  will  shrink  in  the  radial  direction,  which  implies  a  negative  displacement 
Ur.  On  the  other  hand,  the  thermal  expansion  mismatch  between  fiber  and  matrix  will 
induce  a  large  compressive  stress  in  the  fiber,  which  forces  the  fiber  to  expand  in  the 
radial  direction.  For  most  thick  plates,  this  effect  plays  a  much  more  important  role  than 
the  radial  shrinkage  of  the  fiber  itself.  It  is  interesting  to  note  that  in  the  case  of  a  thin  plate 
(large  ratio  of  a/h),  the  displacement  presents  a  negative  value  at  the  center  region  of  the 
plate  (see  Fig.  5.1  for  a/h  =0.3  case).  This  result  seems  reasonable  if  one  more  closely 
examines  the  limit  case.  Suppose  the  plate  is  very  thin  (i.e.,  large  a/h  case).  A  lesser 
compressive  stress  will  be  induced  in  the  fiber  to  expand  it  in  the  radial  direction  and 
the  shrinkage  itself  in  the  radial  direction  due  to  the  decrease  of  the  temperature  will  be  the 
governing  factor,  which  leads  the  displacement  to  be  a  negative  value  over  a  wide  span  of 
the  center  region.  In  the  vicinity  of  the  free  surface,  the  slope  of  each  curve  increases 
rapidly  because  of  the  presence  of  a  stress  singularity  at  z=h,  which  will  be  discussed  later. 
It  is  also  interesting  to  note  that  as  the  a/h  increases  (the  thin  plate  case),  the  curve  around 
the  free  surface  tends  to  smooth.  Because  as  the  plate  tends  to  be  thinner,  the  interaction  of 
the  fiber  and  the  matrix  in  the  axial  direction  tends  to  be  weaker,  which  implies  that  the 
displacement  Ur  tends  to  be  a  constant  along  the  thickness. 

Displacement  w  is  shown  in  Figs.  5.3  and  5.4.  These  curves  start  from  zero  (center 
layer  must  be  zero  because  of  the  symmetry  in  z)  and  they  are  almost  the  linear  functions  of 
z.  Only  in  the  vicinity  of  the  free  surface,  these  curves  present  nonlinear  properties.  As  the 
ratio  of  a/h  decreases  (thick  plate  case),  the  value  of  w  becomes  much  larger  in  comparison 
to  those  obtained  for  large  a/h  ratios.  In  the  case  of  a/h=0.0i,  for  instance,  suppose  the 
inclusion  is  made  of  the  same  material  as  the  matrix  (continuous  plate).  No  mismatch 
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occurs  under  the  change  of  temperature  and  the  value  of  this  displacement  w  at  free 
surface 

z=h  will  be  found  easily  as  -100.  From  Fig  5.3,  however,  the  curve  of  a/h=0.01  attains  the 
value  of  -90.  This  difference  suggests,  therefore,  that  the  inclusion  tries  to  resist  matrix 
shrinkage  along  the  axial  direction. 

The  plot  of  the  stress  a„  in  Fig.  5.5  shows  that  each  curve  looks  parallel  at  a  wide  span 
of  the  center  region  and  there  is  only  a  slight  change  in  the  magnitude  of  the  stress.  Near 
the  surface  of  the  plate,  the  stress  jumps  suddenly  to  a  large  value.  In  view  of  some 
previous  work  (Folias,  1989),  a  stress  singularity  is  present  in  the  neighborhood  of  the 
intersection  of  the  fiber  and  the  free  surface,  i.e.,  at  r=a  and  z=h.  Therefore  the  numerical 
results  shown  in  Fig.  5.5  may  not  be  reliable  near  z=h  because  the  presence  of 
singularities  slows  down  the  convergence  of  the  series  representation  of  the  stresses  near 
the  singularity  point.  Be  that  as  it  may,  the  curves  show  the  trend  of  the  stress  at  z=h  and 
confirm  Folias's  conclusions  [11].  It  may  also  be  noted  that  as  the  shear  moduli  ratio 
changes  from  |i2/M'i=16.67  to  p2/|Ai=6.3  (see  Fig.  5.6)  the  radial  stress  becomes  negative 
over  a  ivide  span  of  the  center  region  and  the  curves  no  longer  appear  to  be  parallel,  but 
mix  together  at  about  four-fifths  the  span  of  the  thickness.  The  reason  for  this  may  be 
explained  as  the  following:  the  radial  stress  o„  induced  by  thermal  contraction  mismatch 
consists  of  two  different  parts.  First,  in  the  radial  direction,  it  is  observed  that  the  inclusion 
presents  a  resistance  to  the  shrinkage  of  the  matrix,  which  leads  to  a  negative  part  of  the 
radial  stress.  On  the  other  hand,  the  inclusion  also  functions  as  a  resistance  to  the 
shrinkage  of  the  matrix  in  the  axial  direction,  which  induces  a  positive  part  of  the  radial 
stress  in  the  matrix.  The  final  sign  of  the  radial  stress  a„  depends  on  which  part  of  the 
effects  is  larger.  From  Figs.  5,5  and  5.6  one  can  see  that  as  the  case  of  shear  moduli  ratio 
increased  or  the  ratio  of  a/h  decreased,  the  second  effect  will  be  the  controlling 
factor  to  the  radial  stress  c„ . 
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The  profiles  of  the  stress  Oeel*)  in  the  matrix  (see  Fig.  5.7)  are  very  close  at  the  center 
region  for  different  ratios  of  a/h.  Then  they  increase  rapidly  near  the  free  surface  due  to  the 
presence  of  a  stress  singularity.  Fig.  5.8  shows  another  case  representing  a  different 
material.  Both  plots  show  that  the  stress  Gee^U  attains  positive  values  for  various  a/h  and 
)i2/p.i  ratios.  This  result  is  expected  because  both  kinds  of  effects  explained  for  the  radial 
stress  Orr  above  are  also  valid  for  the  stress  Oee(*).  The  only  difference  for  the 
circumferential  stress  Gee(i)  is  that  both  parts  of  the  effects  lead  to  positive  values. 

Perhaps  it  may  be  appropriate  here  to  note  that  previous  research  by  Delale  (1988) 
focused  on  a  2-D  consideration  [36].  As  for  the  plane  stress  case,  the  stress  Oee^i)  presents 
a  positive  value  and  the  stress  Gn  presents  the  same  magnitude  but  negative  value  as  the 
stress  G0e(i)  in  his  paper.  Both  the  stresses  are  constants  along  the  thickness  of  the  plate. 
So  the  conclusion  is  clear,  that  is,  a  through  the  thickness  crack  will  initiate  at  the 
fiber/matrix  interface  in  the  radial  direction  as  the  temperature  cools  down  to  a  certain 
degree.  Based  on  his  solution,  if  the  plate  is  thick  enough,  here  a  real  important  fact  has 
been  neglected,  that  is  the  effect  of  mismatch  of  the  fiber  and  matrix  in  the  axial  direction. 
On  the  basis  of  3-D  considerations  as  discussed  in  this  thesis,  a  crack  will  most  likely 
initiate  at  the  stress  singularity  point  near  the  free  surface  and  compare  the  plots  of  stress  G[t 
with  the  plots  of  stress  Gee(0  from  Figs.  5.5  to  5.7,  it  is  difficult  to  conclude  which  stress 
is  the  maximum  because  it  depends  on  the  ratios  of  a/h  and  -  One  may  focus  on  the 
central  region  of  the  plate,  in  the  case  of  a/h=0.01  and  p2/lti=16.67,  the  stress  G„  is  larger 
than  the  stress  Geeli).  In  the  case  of  a/h=0.05  and  the  case  of  a/h=0.3,  however,  the  stress 
Gee^^l  becomes  the  maximum  stress.  This  means  that  the  direction  of  the  crack  initiation  is 
a  function  of  the  ratio  a^  and  of  the  shear  moduli  rallc  I'Vlti- 

The  plots  seen  in  Fig.  5.9  and  Fig.  5.10  show  a  sharp  increase  of  the  stress  Gj^d)  for 
the  plate  in  the  vicinity  of  z=h.  As  the  ratio  of  a/h  decreases,  the  central  region  of  the  plate 
may  not  take  vertical  load  (see  curves  with  a/h=0.01).  As  the  ratio  of  a/h  increases,  the 
stress  level  also  increases.  In  addition,  after  comparing  the  plot  of  Fig.  5.9  with  the  plot  of 
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Fig.  5.10,  one  can  note  that  the  curves  corresponding  to  1=16.67  show  a  higher  stress 
level  than  those  corresponding  to  p.2/|ii=6.3.  The  implication  of  this  is  that  a  higher  ratio 
of|i2^p.i  will  induce  higher  mismatch  and  as  a  result  it  will  translate  a  larger  vertical  stress  to 
the  plate. 

The  shear  stress  Xrz  (see  Fig.  5.1 1-  5.12)  starts  from  zero  at  the  center  of  the  plate  and 
then  increases  in  absolute  value  as  one  travels  along  the  interface.  A  sudden  jump  of  this 
stress  near  the  point  z=h  is  also  observed  and  may  very  well  suggest  why  a  crack  is  most 
likely  to  initiate  at  such  regions.  The  negative  sign  means  that  the  shear  stress  will  always 
be  opposing  the  direction  of  the  axial  stress  Ozz(2)  of  the  inclusion  in  order  to  satisfy  the 
equilibrium  in  the  axial  direction  for  the  fiber. 

The  vertical  stress  in  the  fiber  presents  a  negative  value  in  any  circumstance 
because  the  shrinkage  of  the  matrix  in  the  axial  direction  is  much  larger  than  the  shrinkage 
of  the  fiber  and  it  will  pull  the  fiber  down  in  order  to  satisfy  the  perfect  boundary  condition 
at  the  interface.  The  small  ratio  of  a/h  will  induce  a  large  axial  compressive  stress  in  the 
fiber,  which  increases  further  as  the  ratio  of  Ii2/M^i  also  increases.  The  maximum 
compressive  axial  stress  is  located  at  the  center  of  the  fiber  length  (i.e.,  z=0).  It  is 
interesting  to  compare  the  curves  of  ChW)  with  the  curves  of  0^(2) ,  which  will  show  us  a 
dramatic  situation;  as  the  thickness  of  the  plate  increases  (i.e.,  small  a/h),  the  axial  stress  in 
the  matrix  decreases  while  the  axial  stress  in  the  fiber  increases  rapidly,  which  suggests, 
therefore,  that  the  gap  of  the  stress  discontinuity  at  the  interface  will  increase.  Thus,  the 
ratio  of  a/h  may  be  used  to  control  the  stress  distribution  between  the  fiber  and  the  matrix. 

The  comparisons  of  current  solutions  with  the  solutions  of  Delade  [36]  (based  on  2D 
plane  stress  consideration)  arc  shown  in  Fig.5.15  to  Fig.  5.18.  One  is  able  to  extract 
direedy  from  those  plots  that  as  the  thickness  of  the  plate  decreases  (i.e.,  a^  increases),  the 
2-D  solution  is  recovered  for  both  radial  stress  and  circumferential  stress.  This  conclusion 
can  be  proved  easily  from  the  equations  (4.3)-(4.8).  Readers  may  note  that  as  the  thickness 
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of  plate  increases,  the  ratio  of  radial  stress  a,r(3D)/<Jrr(2D)  decreases  (<1)  whereas  the  ratio 
of  circumferential  stress  O00(3D)/o^(2D)  increases  (>1). 
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5.2  Periodic  Fibers  Model 

In  as  much  as,  in  practice,  the  composite  plate  is  not  embedded  with  only  one  fiber,  it  is 
desirable  now  to  extend  the  model  to  a  more  general  case  where  an  array  of  periodic  fibers 
is  introduced  to  embed  into  the  matrix.  As  for  this  periodic  fibers  model,  all  the 
displacements  and  stresses  will  be  6-dependent  due  to  the  interaction  of  adjacent  fibers.  A 
new  ratio  a/b  (the  ratio  of  fiber  radius  to  separation  distance  between  adjacent  fibers)  is 
introduced  here  to  correspond  to  the  fiber  volume  fraction  Vf .  The  displacement  ue  and  the 
stresses  t,©  ,  Xqz  are  no  longer  zero  here  but  they  are  still  not  the  important  factors 
compared  with  the  other  ones.  Therefore,  they  will  not  be  discussed  in  this  thesis. 

The  shear  moduli  ratio  used  here  is  fixed  at  16.67  throughout  the  discussion. 
Several  different  geometrical  ratios  a/b  are  chosen:  they  are  a/b=0.55,  a/b=0.45,  a/b=0.32 
or  a/b=0.2.  Another  geometrical  ratio  a/h  is  taken  as  before,  that  is,  0.01,  0.05  and  0.3. 
The  plots  are  presented  in  six  stages.  In  the  first  stage,  the  displacement  and  stress  fields 
are  shown  as  the  function  of  0  and  are  compared  with  different  ratios  of  a/b.  In  the  second 
stage,  the  effect  of  different  plate  thickness  is  investigated  where  the  displacement  and 
stress  fields  are  functions  of  z.  In  the  third  part,  the  stresses  vs  a/b  are  examined.  Again, 
the  stress  field  is  presented  as  a  function  of  z  in  the  fourth  part,  but  here  it  is  used  to 
compare  with  different  ratios  of  a/b  instead  of  the  different  ratios  of  a/h.  In  the  fifth  stage, 
the  displacement  and  stress  fields  are  examined  as  a  function  of  the  radial  coordinate  r. 
Finally,  the  effect  of  different  shear  moduli  \i2l\i\  is  investigated. 

Fig.  5.19  shows  the  displacement  Ur  vs  0  at  the  interface  and  along  the  center  of  the 
plate.  The  maximum  deformation  occurs  at  6=0.  It  is  interesting  to  note  that  as  the  ratio  a/b 
increases,  the  value  of  the  displacement  increases  at  0=0  and  decreases  at  0=45.  As  the 
ratio  a/b  decreases  (the  distance  of  adjacent  fibers  increases),  the  curve  tends  to  be  constant 
in  the  0  direction,  which  seems  reasonable  if  one  thinks  of  a  limit  case  a/b->0,  the  result  of 
the  one  fiber  model  is  recovered  (0-independence). 
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Figure  5.19.  Displacement  u,  at  r=a,z=0  vs  0  for  periodic  fibers  model 
where  vi  =0.34, V2=0.22,ji2/p.i=l  6.67  and  a/h=0.05 
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The  displacement  w  at  the  free  surface  is  shown  in  Fig.  5.20.  The  curves  are  close  at 
0=0,  then  separate  in  the  other  side  (0=45).  The  maximum  value  in  magnitude  is  also  at 
0=0  and  as  the  ratio  a/b  decreases,  the  displacement  also  tends  to  be  independent  of  0. 

The  value  of  radial  stress  (s„  (Fig.  5.21)  increases  at  one  side  (0=0)  and  decreases  at  the 
other  side  (0=45)  as  a/b  increases.  The  circunoferendal  stress  Oee(i)  shown  in  Fig.  5.22 
presents  a  similar  condition.  The  maximum  stress  for  both  o„  and  OeeO)  is  located  at  0=0 
for  any  different  ratios  of  a/b.  Fig.5.23  shows  us  the  profiles  of  axial  stress  OzzO)  vs  0. 
All  curves  merge  together  at  about  0=45.  The  maximum  axial  stress  of  OnO)  is  also  located 
at  0=0  but  the  axial  stress  is  small  in  comparison  to  the  radial  stress  0^  and  the 
circumferential  stress  OeeO). 

Figs.  5.24  to  5.34  show  the  displacement  and  stress  fields  vs  ?/h  at  the  interface  of  fiber 
and  matrix  for  both  0=0  and  0=45  cases  where  the  ratio  a/b  is  fixed  at  0.55.  As  it  has  been 
shown  in  Figs.  5.19  to  5.23,  the  maximum  stress  is  located  at  0=0.  Thus  in  the  subsequent 
discussion,  one  need  only  concentrate  on  the  case  0=0. 

The  plots  of  the  displacement  Ur  shown  in  Figs.  5.24  and  5.25  provide  us  different 
profiles  for  the  displacement  at  0=0  and  0=45.  The  variation  of  the  curves  seems  more 
regular  in  the  position  of  0=0  than  the  curves  of  0=45  because  the  curves  of  the  former  one 
cross  at  about  the  same  point  (z/h=0.8)  where  the  value  of  each  curve  is  almost  zero.  As  the 
ratio  a/h  decreases,  the  value  of  the  displacement  increases  at  the  central  region 
(0<z/h<0.8),  and  at  the  region  of  0.8<z/h<l,  it  increases  again  but  in  a  negative  sense.  It 
may  also  be  noted  that  for  a  large  ratio  of  a^  (the  case  of  a/h=0.3),  the  displacement  Ur 
vanishes  in  a  wide  span  from  the  central  region  and  one  is  able  to  extract  that  this 
displacement  tends  to  be  a  negative  value  as  the  a/h  increases  again.  The  displacement  w 
(Figs.  5.26  and  5.27),  on  the  other  hand,  looks  more  regular  than  other  components 
because  it  presents  a  linear  relationship  along  the  z  direction  except  the  values  near  the  free 
surface. 

Figs.  5.28  and  5.29  show  the  curves  of  radial  stress  ct„  at  the  interface  corresponding 
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Figure  5.20.  Displacement  w  at  r=a,  z=h  vs  0  for  periodic  fibers  model 
where  vi  =0.34, v2=0.22,H2/Hi=1  6.67  and  a/h=0.05 
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Figure  5.23.  Stress  diP)  at  r=a,  z=0  vs  0  for  periodic  fibers  model 
where  vi=0.34,V2=0.22,p.2/M^i=16.67  and  a/h=0.05 
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Figure  5.24.  Displacement  Ur  at  T^=a,0=O  vs  z/h  for  periodic  fibers  model 
where  Vi=0.34,v2=0.22,^i2/M.i=16.67  and  a/b=0.55 
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Figure  5.25.  Displacement  u,  at  r=a,6=7C/4  vs  z/h  for  periodic  fibers  model 
where  vi=0.34,v2=0.22,p.2/lti=16.67  and  a/b=0.55 
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Figure  5.26.  Displacement  w  at  r=a,0=O  vs  z/h  for  periodic  fibers  model 
where  vi=0.34,v2=0.22,p2/M-i=16.67  and  a/T)=0.55 


Figure  5.27.  Displacement  w  at  r=a,0=n/4  vs  z/h  for  periodic  fibers  model 
where  vi=0.34,V2=0.22,Li2/iii=16.67  and  a/b=0.55 


Figure  5.28.  Stress  CJ„  at  r=a,6=0  vs  z/h  for  periodic  fibers  model 
where  vi  =0.34, V2=0.22,p2/4i=l  6.67  and  a/b=0.55 


Figure  5.29.  Stress  Oa  at  r=a,0=7C/4  vs  z/h  for  periodic  fibers  model 
where  vi=0.34,v2=0.22,ii2/lti=16-67  and  a/b=0.55 


Figure  5.30.  Stress  099(0  at  r=a,0=O  vs  z/h  for  periodic  fibers  model 
where  vi=0.34,V2=0.22,u.2/Ui=16.67  and  a/b=0.55 


Figure  5.31.  Stress  CqqW  at  r=a,0=7C/4  vs  z/h  for  periodic  fibers  model 


where  vi=0.34,V2=0.22,Li2/Hi=16.67  and  a/b=0.55 


Figure  5.32.  Stress  CuW  at  r=a,0=O  vs  z/h  for  periodic  fibers  model 
where  vi  =0.34, V2=0.22,^2/4i=l  6.67  and  a/b=0.55 


Figure  5.33.  Stress  Xrz  at  r=a,8=0  vs  z/h  for  periodic  fibers  model 
where  vi=0.34,V2=0.22,ii2/Ui=16.67  and  a/b=0.55 


Figure  5.34.  Stress  Czz®  at  r=a,0=O  vs  z/h  for  periodic  fibers  model 
where  vi*0.34,V2=0.22,U2/Ui=l6.67  and  a/b=0.55 
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to  two  different  locations  0=0  and  0=45.  The  radial  stress  increases  with  the  decrease  of 
a/h  in  both  directions.  Readers  can  note  from  Fig.  5.28  that  as  a/h  reaches  a  certain  value 
(small  a/h),  all  the  curves  tend  to  merge  together  at  the  central  region  of  the  plate,  which 
reveals  that  as  the  thickness  of  the  plate  increases  to  a  certain  degree,  the  radial  stress  will 
stop  increasing  and  reach  a  limit  value  in  the  central  region  of  the  plate. 

The  circumferential  tress  aee(i)  is  shown  in  Figs.  5.30  and  5.31,  where  there  is  only  a 
slight  change  of  the  magnitude  at  the  central  region.  However,  the  circumferential  stress 
jumps  suddenly  to  a  large  value  near  the  plate  surface.  A  special  case  is  the  curve  of 
a/h=0.01  at  0=0.  In  this  case,  the  circumferential  stress  decreases  suddenly  from  the  point 
of  about  z/h=0.7  then  returns  back  rapidly  from  about  z/h=0.95  to  a  large  value.  The  axial 
stress  Ozz(i)  (Fig.5.32)  for  the  case  of  a/h=0.0l  has  a  phenomenon  similar  to  OeeO) .  The 
shear  stress  Trz  (see  Fig.  5.33)  increases  in  magnitude  almost  linearly  at  the  central  region 
for  every  ratio  a/h  then  increases  fast  near  the  free  surface  where  the  stress  singularity  is 
located.  Finally,  the  vertical  stress  of  the  inclusion  azz(2)  is  examined  in  Fig.  5.34. 
Comparing  the  plot  with  the  plot  in  Fig,  5.13,  one  can  note  that  both  plots  look  similar  , 
which  reveals  that  the  effect  of  the  ratio  a/b  to  the  stress  <Tz2(2)  is  negligible. 

Figs.  5.35  to  5.37  show  the  variation  of  the  stress  profiles  a„ ,  099(1)  and  Ozz(')  vs 
a/b  for  the  case  of  a/h=0.01,  respectively.  Each  of  these  plots  contains  two  curves  that 
represent  two  different  points  at  0=0  and  0=45  in  the  center  of  the  interface.  In  Fig.  5.35, 
the  radial  stress  Orr  decreases  in  both  directions  as  the  ratio  a/b  increases.  On  the  other 
hand,  the  curves  of  the  circumferential  stress  shown  in  Fig.  5.36  are  separated  in  different 
directions  with  the  increase  of  a^.  It  is  of  interst  to  note  that  as  the  distance  of  the  adjacent 
fibers  decreases,  parts  of  the  stress  099(1)  will  transfer  from  0=45  to  0=0  to  add  the 
magnitude  of  the  circumferential  stress  099(1)  at  0=0.  The  axial  stress  Ozz(i)  vs  a/b  is 
attained  in  Fig.5.37.  It  is  observed  that  this  axial  stress  increases  fast  at  0=0  with  the 
increase  of  ratio  a/b  whereas  it  changes  only  slightly  at  0=45. 

Figs.  5.38  to  5.42  show  the  profiles  of  those  normal  stresses  a„ ,  099(1)  and  Ozz(l) 


Figure  5.38.  Stress  o„  at  r=a,  0=0  vs  z/h  for  periodic  fibers  model 
where  vi=0.34,V2=0.22,p.2/lti=16.67  and  a/h=0.01 


Figure  5.39.  Stress  Orr  at  r=a,  0=7t/4  vs  z/h  for  periodic  fibers  model 
where  vi=0.34,V2=0.22,^2/Ui=16.67  and  a/h=0.01 


Figure  5.40.  Stress  at  r=a,  0=0  vs  z/h  for  periodic  fibers  model 
where  vi  =0.34, V2=0.22,U2/Ui=l  6.67  and  a/h=0.01 


Figvire  5.41.  Stress  Oee(i)  at  r=a,  0=7t/4  vs  z/h  for  periodic  fibers  model 


where  vi=0.34,V2=0.22,H2/4i=16.67  and  a/h=0.01 


Figure  5.42.  Stress  OzzCO  at  r=a,  0=0  vs  z/h  for  periodic  fibers  model 
where  vi=0.34,V2=0.22,fi2^1ii=16.67  and  a/h=0.01 


97 


again  but  using  different  ratios  of  a/b  instead  of  different  ratios  of  a/h  in  order  to  observe 
more  clearly  the  effect  of  the  ratio  a/b  on  the  stresses  along  the  thickness  and  at  the 
interface.  In  the  position  of  0=0,  the  effect  of  ratio  a/b  is  small  to  the  radial  stress  in  the 
central  region  (see  Fig.  5.38).  In  the  position  of  0=45,  however,  the  different  ratios  of  a/b 
separate  the  radial  stress  as  a  group  of  parallel  curves  (see  Fig.  5.39).  The  same 
phenomenon  can  also  be  seen  in  the  profile  of  the  circumferential  stress  (Figs.  5.40  and 
5.41).  On  the  other  hand,  the  axial  stress  with  different  ratios  of  a/b  at  0=0  is  shown 
in  Fig.  5.42.  The  stress  is  redistributed  by  different  ratios  a/b  along  the  thickness.  As  a/b 
increases,  the  stress  in  the  central  region  increases  but  the  stress  near  the  free  surface 
decreases  to  a  negative  value,  then  goes  back  again  to  a  big  value. 

So  far,  all  of  the  plots  discussed  above  are  focused  on  the  interface  of  the  fiber  and  the 
matrix  (r=a)  where  a  maximum  stress  is  presented.  The  displacement  and  stress  distribution 
along  the  radial  direction  is  also  important  for  us  in  order  to  understand  the  material 
properties.  For  this  reason,  the  plots  (Figs.5.43  to  5.51)  are  made  here  to  examine  how 
displacements  and  stresses  will  decay  as  a  function  of  the  radial  distance  .  All  of  the  curves 
discussed  here  are  within  the  region  begining  at  the  interface  r=a  and  ending  at  the 
symmetric  line  B-C  (see  Fig.  2.2). 

The  displacement  w  at  the  free  surface  z=h  is  shown  in  Figs.  5.43  and  5.44  as  a 
function  of  the  radial  coordinate  r.  At  0=0,  the  displacement  looks  more  complicated  than 
the  displacement  at  0=45  because  the  curve  at  0=0  changes  from  a  concave  shape  to  a 
convex  shape  as  the  ratio  a/b  increases  and  the  curves  at  0=45  look  more  smooth  and 
regular.  The  values  of  the  displacement  at  the  interface  r=a,  which  is  smaller  than  100  in 
magnitude,  reveal  that  the  inclusion  resists  the  shrinkage  of  the  matrix. 

The  radial  stress  profiles  at  the  center  of  the  plate  (z=0)  are  shown  in  Figs.  5.45  and 
5.46  along  the  two  radial  directions  of  0=0  and  0=45,  respectively,  for  different  ratios  of 
a/b.  The  stress  tends  to  decrease  as  the  distance  away  from  the  interface  increases.  It  may 
also  be  noted  that  the  stress  along  the  direction  0=0  decays  faster  than  the  stress  in  the  other 
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Figure  5.44.  Displacement  w  at  z=h,  e=ji/4  vs  (r-a)/(0.707b-a)  for  periodic  fibers  model 
where  vi=0.34,V2=0.22,H2/|Ai=16.67  and  a/h=0.01 


Figure  5.45.  Stress  a„  at  z=0, 0=0  vs  (r-a)/(b-a)  for  periodic  fibers  model 
where  vi=0.34, ¥2=0.22,^2/^1=16.67  and  a/h=0.01 
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Figure  5.46.  Stress  a„  at  z=0, 0=Ji/4  vs  (r-a)/(0.707b-a)  for  periodic  fibers  model 
where  vi=0.34,V2=0.22,42/lti=16.67  and  a/h=0.01 
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Figure  5.47.Stress  099(0  at  z=0, 0=0  vs  (r-a)/(b-a)  for  periodic  fibers  model 


where  vi=0.34,V2=0.22,H2/lii=l6.67  and  a/h=0.01 
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Figure  5.48.Stress  aged)  at  z=0, 0=7C/4  vs  (r-a)/(0.707b-a)  for  periodic  fibers  model 

where  vi=0.34,V2=0.22,p.2/M.i=16.67  and  a/h=0.01 


Figure  5.49.Stress  Ozz^i)  at  z=0,  0=0  vs  (r-a)/(b-a)  for  periodic  fibers  model 
where  vi=0.34,v2=0.22,p2/lii=16.67  and  a/h=0.01 
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Figure  5.51. Stress  at  z=0,  vs  r/a  for  periodic  fibers  m 

where  vj  =0.34, V2=0.22,|X2/lti=l 6.67  and  a/h=0.01,a/b=0.f 
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direction  (0=45).  A  similar  trend  also  prevails  for  the  circumferential  stress  aeeO)  (see 
Figs.  5.47  and  5.48).  As  for  the  stress  OzjCO  (see  Figs.  5.49  and  5.50),  on  the  other 
hand,  one  can  note  that  both  groups  of  curves  along  different  radial  directions  (i.e.,  0=0 
and  0=45)  approach  to  the  value  of  zero . 

Fig.  5.51  shows  the  axial  stress  in  the  fibers  vs  r/a.  At  the  interface  r=a,  the 
magnitude  of  the  axial  stress  increases  along  0=45  and  decreases  along  0=0.  It  seems 
reasonable  if  one  compares  these  curves  with  the  curve  of  a/b=0.55  in  Fig.  5.20.  The 
shrinkage  in  vertical  direction  at  the  interface  is  resisted  much  more  at  0=45  than  at  0=0  (if 
without  the  resistance,  the  value  of  the  shrinkage  is  -100),  which  induces  a  higher 
compressive  stress  Czz(2)  at  the  point  0=45  than  at  the  point  of  0=0. 

Finally,  the  stresses  <s„  and  CTeeCO,  which  are  two  of  the  most  important  stresses  under 
consideration,  are  shown  again  as  a  function  of  shear  moduli  ratio  112/^1  in  Figs.  5.52  and 
5.53,  respectively.  Starting  from  =0.01,  one  can  travel  along  the  curves  and  note 
that  both  radial  stress  and  circumferential  stress  are  close  to  zero  in  the  region  of 
small  shear  moduli  ratios.  This  result  seems  reasonable  if  one  examines  the  limit  case  of  a 
plate  with  periodic  holes.  From  Fig.5.52,  the  reader  may  note  that  the  radial  stress 
appears  to  be  low  for  a  wide  span  of  the  ratio  Beyond  the  ratio  of  =10,  it 

increases  rapidly  .  From  those  four  curves  with  different  parameters,  one  is  able  to  extract 
directly  that  the  effect  of  the  ratio  a/h  to  the  radial  stress  is  much  higher  than  the  effect  of  the 
ratio  a/b.  The  circumferential  stress  OqqW,  on  the  other  hand,  has  a  large  slope  in  the 
region  for  small  ratios  of  112/^1 ,  while  the  rate  of  increase  slows  down  in  the  region  of 
high  ratios  of  H2/^i-  A  reverse  condition  appears  for  this  circumferential  stress  (see  Fig. 
5.53),  that  is,  the  effect  of  the  ratio  a/b  to  this  circumferential  stress  is  higher  than  that  of 
the  ratio  a/h.  Thus,  the  following  conclusion  may  be  drawn  that  in  the  case  of  small  ratios 
of  ,  a  radial  crack  is  most  likely  to  initiate  first,  while  in  the  case  of  large  ratios  of 
IJ.2/II1  a  fiber-matrix  interface  crack  is  most  most  likely  to  develop. 


CHAPTER  6 


CONCLUSIONS  AND  RECOMMENDATIONS 

6.1  Conclusions 

In  view  of  the  numerical  results  discussed  above,  the  following  conclusions  may  be 
made: 

1.  The  thickness  and  distance  of  adjacent  fibers  of  a  plate  ,  as  well  as  the  material 
properties,  play  a  fundamental  role  in  the  failure  mechanism  of  a  plate  with  one  or  periodic 
cylindrical  inclusions.  More  specifically,  the  geometric  parameters  a/h,  a/b  and  the  shear 
moduli  ratio  ^  three  important  ratios  that  not  only  control  the  location  where  a 
failure  initiation  may  occur,  but  may  also  determine  the  direction  where  the  crack 
propagates. 

2.  The  value  of  shear  moduli  ratio  42/^1  may  also  be  used  to  predict  the  direction  of  a 
crack  initiation.  As  for  a  small  ratio  of  |X2/lti  (<10),  a  radial  crack  may  initiate  first  whereas 
for  a  large  ratio  of  112^111  (>I3),  a  crack  along  the  fiber-matrix  interface  may  initiate  first 
then. 

3.  The  ratio  of  thermal  expansion  coefficient  a(2)/a(i)  presents  a  linear  action  to  the 
displacen^nt  and  stress  fields. 

4.  Under  the  cooling  process  as  discussed  in  this  thesis,  the  failure  (debonding, 
slippage  and  breakage)  will  most  likely  be  induced  by  the  effect  of  the  radial  stress  On  or 
the  circumferential  stress  Oee^^)  •  One  can  predict,  however,  a  conclusion  based  on  the 
information  discussed  in  Chapter  5,  if  the  plate  is  subjected  to  a  high  temperature.  In  the 
case  of  an  airplane  with  a  high  speed,  a  crack  may  initiate  at  the  central  region  of  the  fiber 
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because  in  this  circumstance,  the  axial  stress  azz(2)  in  the  fiber  tends  to  be  a  large  tensile 
stress  that  will  control  the  failure  of  the  plate. 

5.  The  comparison  of  the  current  result  with  the  result  of  Delale  (1988)  indicates  that  the 
current  solution  can  predict  the  interfacial  shear  stress  quite  satisfactorily  throughout  the 
interface,  whereas  the  2-D  solution  obtained  by  Delale  [36]  neglects  the  effect  of  the 
mismatch  along  the  axial  direction,  which  has  been  proved  to  be  a  very  important  factor  to 
the  stress  field  for  a  thick  plate. 

6.  An  edge  effect  is  confirmed  to  exist  in  the  vicinity  of  the  intersection  of  the  interface 
(r=a)  and  the  free  surface  (z=h)  where  the  presence  of  a  stress  singularity  may  result  in 
crack  initiation. 

7.  Finally,  to  avoid  unnecessary  damage  induced  by  the  residual  stress,  an  optimal 
design  of  fiber-reinforced  composite  material  may  be  recommended  to  manufacturing 
engineers; 

(a) .  The  shear  moduli  ratio  (fiber  to  matrix)  is  a  factor  that  needs  to  be 
optimized.  A  high  value  of  this  ratio  may  cause  a  fiber/matrix  interface  crack  (see  Fig. 
5.52)  and  therefore  should  be  avoided  in  most  circumstances. 

(b) .  The  fiber  volume  fraction  Vf  (i.e.,  the  ratio  of  a/b)  as  well  as  the  ratio  a/h  are 
also  factors  that  need  to  be  optimized.  The  maximum  stress  may  change  between  the  radial 
stress  and  the  circumferential  stress  with  the  change  of  the  ratios  a/b  and  a/h.  Thus,  two 
types  of  failure  immediately  come  to  mind:  a  fiber/matrix  interface  crack  and  a  radial  matrix 
crack.  One  is  able  to  deduce  that  the  minimum  critical  residual  stress  occurs  at  the  position 
where  the  radial  stress  g„  equals  the  circumferential  stress  OeeO)  •  The  ratios  of  a/b  and 
a/h  may  then  be  thought  of  as  the  optimal  ratios  to  reduce  the  residual  stress  effect. 

(c) .  The  thermal  expansion  coefficients  ratio  a(2)/a(i)  presents  a  linear  action  to  the 
residual  stresses.  It  is  optimal  to  use  the  materials  of  fiber  and  matrix  with  close  thermal 
expansion  coefficients  to  minimize  the  effect  of  the  residual  stresses. 
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6.2  Recommendarions  for  Further  Work 

The  present  analysis  does  give  some  important  information  that  can  be  used  in  the  study 
of  composite  materials.  However,  much  work  still  remains  to  be  done  in  the  future. 
General  speaking,  Folias'  3-D  solution  can  be  developed  to  solve  more  complicated 
problems  if  one  knows  how  to  modify  the  program  to  remedy  the  ill  conditions.  Based  on 
the  author's  experience  and  knowledge,  the  following  problems  need  to  be  addressed  in  the 
future: 

(1) .  Perfect  bonding  and  isotropic  material  of  both  fiber  and  matrix  were  assumed  to 
prevail  in  this  paper.  However,  a  large  mismatch  has  been  found  in  the  interface  along  the 
thickness.  A  slip  in  the  axial  direction  may  occur.  So  a  new  model  that  relaxes  the 
assumption  of  perfect  bonding  to  a  partial  debonding  case  must  be  dealt  with.  On  the  other 
hand,  anisotropic  properties  exist  in  some  materials,  such  as  the  carbon  fiber.  So  the 
solution  at  hand  may  be  developed  to  a  more  general  case  that  is  valid  for  an  anisotropic 
material. 

(2) .  Some  pre-existing  cracks  can  be  introduced  to  the  investigation.  The  particular 
solution  in  the  method  needs  to  be  modified  to  handle  the  mixed  condition  at  crack  region. 

(3).  All  of  the  previous  research  is  limited  to  the  "  micro-approach,"  which  considers  the 
micromechanical  behavior  of  only  one  ply  of  laminate.  In  practice,  the  interaction  of  the 
laminates  and  edge  effect  between  the  plies  with  different  phases  is  also  very  important. 
How  to  develop  a  model  that  can  be  used  to  solve  the  multiply  composites  is  an  important 
consideration  for  future  research. 


APPENDIX  A 

THE  SMALL  CONTRIBUTIONS  OF  THE  MODIFIED  BESSEL 

FUNCTION  I  TO  THE  STRESS  FIELD  IN  THE  MATRIX 
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In  this  thesis,  only  the  Bessel  function  K  is  used  to  match  the  stress  field  in  the  matrix  (r 
>  a  region).  Because  the  residual  stresses  are  induced  by  the  mismatch  of  fiber/matrix 
interface  (r=a),  the  stresses  in  the  matrix  decrease  along  the  radial  direction  r,  which 
matches  the  trend  of  the  modified  Bessel  function  K.  The  other  modified  Bessel  function, 
I,  however,  goes  to  infinity  as  r  increases;  therefore  it  is  expected  that  the  coefficient  that 
multiplies  this  function  must  decrease  rather  rapidly. 

In  order  to  confirm  that  the  effect  of  adding  the  term  of  the  modified  Bessel  function  I  to 
the  matrix  can  be  neglected,  an  efficient  procedure  is  as  follows.  Consider  the  linear  system 
of  both  the  governing  equations  and  boundary  conditions,  which  may  be  expressed  as 
follows: 


L[U]-  ) 

(governing  equations) 

(A-1) 

M[U]=0 

(boundary  conditions) 

(A-2) 

U  =  u  +  u*  , 

(A-3) 

where  u  represents  the  part  that  is  used  in  this  thesis  and  u*  represents  the  solution  with 
the  modified  Bessel  function  I.  Then  from  equations  (D-l)-(D-3),  one  gets: 

L  [u+u*  1  =  L  [u]  +  L[u*  ]  =  0  +  L[u*  ]  =  L  [u*  ]  =  0  (A-4) 

M  [u+u*  ]  =  M  [u]  +  M  [u*  ]  =  e  +  M  [u*  ]  =  0  (A-5) 

where  e  is  the  error  of  the  boundary  conditions  obtained  before.  It  may  be  noted  that 
equation  (D-4)  is  satisfied  automatically  and  that  the  only  problem  remaining  is  to  solve 
equation  (D-5),  which  may  be  expressed  as  a  new  linear  system: 


A*  X  =  C  =  -  e 


(A-6) 
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Solving  for  the  constants  x  and  computing  the  additional  stress  field  a* ,  one  can  see 
clearly  that  the  contribution  of  the  modified  Bessel  function  I  to  the  stress  field  in  the  matrix 
is  negligible. 

Figs.D.l  to  D.4  show  us  the  profiles  of  the  ratio  of  this  additional  stress  to  the  stresses 
obtained  before.  Both  the  radial  stress  a„  and  the  circumferential  stress  099  are  shown  as 

a  function  of  r/a  or  z/h,  respectively.  Two  different  sets  of  roots  are  given  in  the  plots  to 
show  the  difference  of  the  value.  It  is  noted  from  these  plots  that  the  contribution  due  to  the 
extra  term  is  less  than  0.5%  of  the  stress  field,  which  was  obtained  previously.  This 
confirms  the  assumption  that  the  contributions  due  to  the  modified  Bessel  function  I  may  be 
neglected  and  still  have  accuracy  up  to  1%.  By  utilizing,  however,  both  solutions,  one 
may  calculate  the  stress  field  to  within  an  error  of  less  than  0.1%. 


(r-a)/(b-a) 


Figure  D.l.  The  ratio  of  stress  at  z=0, 0=7C/4  vs  (r-a)/(b-a)  for  periodic  fibers  model 
where  vi=0.34,  V2=0.22,p2/Hi=16.67,  a/b=0.55  and  a/h=0.05 


Figure  D.2.  The  ratio  of  stress  Oee  at  z=0,  0=7t/4  vs  (r-a)/(b-a)  for  periodic  fibers  model 
where  vi=0.34,  V2=0.22,U2/u.i=16.67,  a/b=0.55  and  a/h=0.05 
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ABSTRACT 


This  paper  presents  several  topics  related  to  the  3D  mechanics  study  of  a  single  fiber 
embedded  into  a  matrix  and  subjected  to  axial  loadings.  Both  the  perfect  bonding  model 
and  the  imperfect  bonding  model  are  formulated  and  numerically  solved  by  a  method 
developed  by  Folias  (1975, 1990).  The  numerical  results  show  that  the  geometric  ratio  a/h 
(fiber  radius  to  plate  half  thickness)  and  shear  moduli  ratio  (Hber  to  matrix)  have 
great  effect  on  the  stress  field.  The  failure  n»chanisms  corresponding  to  these  different 
models  are  found  in  the  paper.  In  addition,  comparisons  with  the  experimental  results  are 
also  made. 
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INTRODUCTION 

The  problem  of  determining  the  stress  distribution  induced  by  inclusions  under  load  has 
been  a  hot  topic  for  about  half  a  century.  Two-dimensional  solutions  (plane  stress  and 
plane  strain)  for  plates  with  perfectly  bonded  circular  inclusions  can  be  found  in  the  papers 
of  Sendeckyj  (1970)  and  of  Yu  and  Sendeckyj  (1974).  A  general  representation  of  the 
solution  of  the  elastic  curvilinear  inclusion  problem  is  presented  by  Sendeckyj.  As  an 
example,  the  elliptical  shape  of  inclusion  is  discussed.  The  discussion  is  limited  to  the  case 
of  an  infinite  matrix.  Later,  the  problem  of  an  unbounded  elastic  matrix  containing  any 
number  of  elastic  inclusions  is  solved  by  Yu  et  al..  Both  of  these  use  the  complex  variable 
formulation.  A  more  practical  model  for  the  mechanical  behavior  of  unidirectional  tiber- 
reinforced  materials  subjected  to  axial  loading  is  examined  by  Bloom  (1967),  where  a 
hexagonal  array  of  perfectly  bonded  filaments  is  assumed.  While  these  qrproaches  are  well 
understood  in  two  dimensions,  they  are  less  so  in  three.  Since  these  solutions  are  based  on 
two-dimensional  considerations,  the  effect  of  the  thickness  of  the  plate  on  the  stress 
distributions  could  not  be  examined. 

Three  dimensional  solutions  to  similar  problems  are  not  fully  investigated  due  to  the 
mathematical  difficulties  involved.  Muki  and  Sternberg  in  ^970  investigated  the  diffusion 
of  an  axial  load  from  a  bar  of  arbitrary  uniform  cross-section  that  is  immersed  in,  up  to  a 
finite  depth,  and  bounded  to  a  semi-infinite  solid  of  distinct  elastic  properties.  Their 
approximate  method  requires  the  radius  of  the  rod  to  be  small  in  comparison  to  its  length. 
Luk  and  Keer  in  1979  investigated  a  very  similar  problem.  The  rod  bar  at  this  time  is 
assumed  to  be  rigid.  Many  plots  of  the  stresses  based  on  numerical  calculation  are  given. 
The  authors  have  also  examined  what  effect  different  parameters  of  die  problem  have  on  the 
stress  field.  It  is  important  to  note  that  all  of  the  above  works  deal  with  one  perfect  isolated 
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inclusion.  On  the  other  hand,  in  a  fiber-reinforced  composite  thousands  of  fibers  can  be 
used  to  construct  one  layer  of  the  laminate.  Similar  to  the  geometry  of  the  problem 
examined  by  Bloom  (1967),  a  three-dimensional  solution  is  achieved  by  using  Boussinesq- 
Papkovitch  potentials  by  Haener  (1967).  However,  only  few  numerical  results  are 
presented  in  the  paper.  Perhaps  the  reascm  for  this  is  the  numerical  complexity  encountered 
because  of  the  double  summadon.  Folias  in  1975  examined  a  method  for  constructing 
solutions  to  some  three-dimensional  mixed  boundary-value  problems  and  applied  it  to  the 
problem  of  a  uniform  extension  of  an  infinite  plate  containing  a  through  the  thickness  line 
crack.  Subsequently  the  general  solution  is  used  to  investigate  some  related  problems.  In 
Penado  and  Folias  (1989),  the  stress  Held  around  a  cylindrical  inclusion  in  a  plate  of 
arbitrary  thickness  is  investigated,  where  a  uniform  tension  is  applied  in  the  plane  of  the 
plate  at  points  far  remote  from  the  inclusion.  Since  the  thickness  of  the  plate  is  no  longer 
assumed  to  be  infinite  or  semi-infinite,  the  results  allow  the  examination  of  veiy  thick  and 
very  thin  plates  and  bridge  the  gap  in  betwwn. 

While  all  above  discussions  are  restricted  to  the  case  of  a  perfect  bonding  condition  at  the 
interface,  there  are  few  analytical  models  which  deal  with  the  imperfect  bonding  problem. 
The  fint  model  assumes  that  the  inclusion  and  the  plate  are  connected  by  an  elastic  spring  at 
the  interface.  Tractions  and  displacements  u  and  v  are  continuous  across  the  interface,  but 
the  vertical  displacement  w  may  not  be  continuous  and  the  difference  Aw  is  assumed  to  be 
proportional  to  the  shear  stress  Xn  at  the  interface.  Papers  by  Lawrence  (1972)  and  Banbaji 
(1988)  approach  the  problem  based  on  shear  lag  analysis  and  the  results  are,  therefore, 
approximate.  Most  recent  work  done  by  Steif  and  Hoysan  (1986)  iq>proaches  the  problem 
based  on  2D  considerations  and  the  case  in  which  the  fiber  and  the  matrix  have  identical 
elastic  properties  is  solved.  Numerical  solutions  for  cases  in  which  die  fiber  and  die  matrix 
have  different  elastic  properties  are  obtained  by  a  finite  element  method.  The  second  model 
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developed  by  Dollar  and  Steif  (1988)  assumes  that  the  transfer  of  load  at  the  interface  is 
described  by  Coulomb  friction.  The  mathematical  models  for  stick,  slip  and  separation 
conditions  are  given  and  the  residual  stress  is  introduced  into  the  discussion.  These  two 
models  have  their  deficiency  which  will  be  seen  later.  Haritos  and  Keer  (1985)  considered 
the  problem  of  a  finite,  rigid  insert  partially  embedded  in  and  adhesively  bonded  to  an 
elastic  half  plane.  The  problem  is  then  formulated  in  terms  of  a  singular  integral  equation 
which  is  solved  numerically.  The  fiber/matrix  debonding  problem  is  solved  by  Gao  et  al. 
(1988)  based  on  the  fracture  energy  balance.  The  friction  at  the  debonded  interface  and  the 
Poisson  contraction  of  the  fiber  are  included.  Interfacial  friction  is  shown  to  have  a 
significant  effect  on  the  debonded  load.  A  similar  work  is  also  carried  out  later  by  Penn  and 
Lee  (1989).  Shear  lag  models  are  used  both  in  these  two  papers  to  calculate  the  stress  and 
the  displacement  fields. 

A  boundary  layer  effect  is  expected  to  prevail  in  the  neighborhood  of  the  intersection  of  the 
interface  and  of  the  plate  surface.  This  matter  has  recently  been  investigated  by  Folias 
(1989)  based  on  3D  consideration,  where  by  assuming  a  certain  form  of  the  solution  he 
was  able  to  use  an  asymptotic  expansion  developed  by  Williams  (1959)  and  to  extract  the 
explicit  displacement  and  stress  fields.  The  stress  field  is  shown  to  be  proportional  to  p-" 
where  a  is  a  weak  singularity  which  depends  on  the  material  properties  and  the  angle  of 
intersection  of  the  fiber  with  the  free  surface. 

On  the  other  hand,  the  experiments  focusing  on  the  interfacial  strengdi  have  also  been  a  hot 
topic  for  quite  a  while.  Experiments  on  a  two  dimensional  model  of  an  aluminium  alloy 
fiber  in  an  Araldite  resin  are  carried  out  by  Tyson  and  Davies  (1965),  the  interfacial  shear 
stress  is  obtained  by  photoelastic  method.  Experiments  conducted  by  Chua  and  Piggott 
(1985)  have  measured  several  important  parameters,  such  as  interfacial  yield  strength  and 
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an  interfacial  work  of  fracture.  The  stress  transfer  and  ftacture  in  single  fiber/epoxy 
composites  have  been  investigated  by  DiBenedetto  et  al.  (1986)  and  Bascom  et  al.  (1986). 
Their  main  efforts  is  to  find  the  cumulative  distribution  of  critical  fiber  lengths,  which  can 
be  used  to  calculate  the  interfacial  shear  strength. 

The  purpose  of  this  analysis  is  to  consider  a  similar  problem  as  that  of  Penado  et  al. 
(1989),  but  the  plate  is  now  subjected  to  a  uniform  tensile  load  along  the  axis  of  the 
cylindrical  inclusion.  As  far  as  the  authors  know,  no  analytical  work  for  the  problem  under 
consideration  has  ever  been  carried  out  The  method  which  will  be  used  in  this  analysis  is 
the  same  as  that  of  Folias  (1975, 1990).  The  problem  will  subsequently  be  extended  to  two 
other  cases:  first  the  plate  surface  will  also  be  allowed  to  carry  a  portion  of  the  load  and 
second  the  interface  will  be  assumed  to  be  imperfectly  bonded.  A  fracture  mechanics 
approach  is  beyond  the  scope  of  current  analysis.  It  is  expected,  however,  that  the  current 
analysis  will  lead  to  the  calculation  of  strain  energy  more  accurately  in  the  future.  In  order 
to  simplify  the  mathematical  complexities  of  our  problem,  the  following  assumptions  will 
be  made: 

(1) .  Both  the  plate  and  the  inclusion  are  made  of  isotropic,  homogeneous  and  linearly 
elastic  materials. 

(2) .  As  a  first  step,  only  erne  isolated  inclusion  is  assumed  to  be  embedded  into  the  plate. 
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EQRMULAUQN  QF  THE  PROBLEMS 

Perfect  Bonding  model  1 

we  consido’  a  body  which  occupies  the  space  Ixkoo,  lykoo,  Izl^  and  contains  two  regions 
with  different  elastic  properties.  The  inclusion  is  a  cylinder  of  radius  r=a  through  the 
thickness  (see  Fig.l).  For  convenience,  the  region  r>a  (  plate  )  and  the  region  r<a 
(inclusion)  are  denoted  by  superscripts  (1)  and  (2),  respectively-  Uniformly  distributed 
tensile  stress  Go  is  applied  both  on  the  bottom  and  the  top  of  the  inclusion  and  all  other 
boundaries  are  assumed  free  of  stress. 

In  the  absence  of  body  forces,  the  governing  equation  (Navier's  equation)  for  the 
displacements  uW,  v(‘),  w<‘)  (i=l,2)  is 


where  is  the  3D  Laplacian  operator,  mi=lA'i,  Vj  is  Poisson’s  ratio  and 


au^‘^ 

=  k=l,2,3 


(2.2) 


The  stress-displacement  reladons  are  given  by  Hooke's  law  as: 

o2=2Hi{  i^ePSik +e®l:  W=1.23  (2.3) 

where  m  are  the  respective  shear  moduli. 


As  to  the  boundary  conditions  one  must  require  that 
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as  r-^oo; 


'^xx^jy^zE 

(2.4) 

(2.5) 

at  z=lhl; 

T'iWi)=<^>=0 

(2.6) 

(2.7) 

II 

(2.8) 

atr=a; 

(2.9) 

(2.10) 

also  the  following  continuity  condition  must  hold 

at  r=0 ;  all  stresses  and  displacements  for 

the  inclusion  must  be  bounded 

(2.11) 

Our  approach  is  to  seek  a  solution  in  the  form: 

u<>)=u<pXi)  +  u<cX») 

(2.12) 

v(»)=sv(PXi)  +  v(cXi) 

(2.13) 

Vy(>)=Vv(pX>)  +  w(cXi) 

(2.14) 

where  the  component  with  the  superscript  (p)  represents  the  particular  solution,  and  the 
comptxient  with  the  superscript  (c)  represents  the  complementary  solution. 
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In  general,  the  particular  solution  is  relatively  easy  to  obtain.  It  must  satisfy  the  governing 
equation  (2.1)  as  well  as  the  boundary  condition  far  away  from  the  inclusion.  For  the 
problem  under  consideration,  the  pardcular  soluticm  in  cylindrical  coordinates  is: 

(1)  for  the  plate; 


(2.15) 

(2.16) 

(2.17) 

(2)  for  the  inclusion: 

(2.18) 

uSf>=0 

(2.19) 

^  ^2\L2  1-V2  1-V2  " 

(2.20) 

(2.21) 

(2.22) 

-(pX2)_..(pX2)  _(pX2)  « 

^10 

(2.23) 

where  Ci  is  a  constant  to  be  determined  later  from  the  boundary  conditions  of  the 
complementary  problem. 

In  view  of  the  particular  solution,  one  needs  to  find  six  complementary  displacements 
v(*),  w(‘)  (i=l,  2),  such  that  they  satisfy  both  the  governing  equation  and  the  boundary 
conditions: 
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at  lzl=h; 

■CW^‘W‘W’^0 

(2.24) 

atr=a; 

o(cXl)  JcX2>^Xl),^)0) 

IT  IT  V  IT 

(2.25) 

.j(cXl)  Jc)(2)  (pXD  (pX2) 

(2.26) 

^cX1),^c)(2)^JX1).^(p)(2) 

(2.27) 

ufrXl).y(cX2)^y(pXl)^<pX2) 

at  ff  ff. 

(2.28) 

u(cX1)„(cX2)_(pX1)^,(pX2) 

“eo  “ee  “^“ee  ^ee 

(2.29) 

„(cX1).„(cX2)_(pX1)^,(pX2) 

(2.30) 

It  may  be  appropriate  at  this  stage  to  point  out  that  a  uniform  tensile  stress  gq  assumed  in 
Fig.l  has  neglected  the  stress  concentration  effect  existing  at  die  edge  region.  Such  a  stress 
concentration  effect  could  be  taken  into  account  if  one  considers  a  non-uniform  tensile 
stress  acting  on  the  fiber  surface  as  shown  in  Fig.  2.  It  will  be  shown  later  that  the  error 
brought  on  by  such  a  simplificadon  is  very  localized  (about  one  fiber  diameter  from  the 
fiber  end).  A  follow-up  paper  on  this  irmter  is  under  preparation. 

Perfect  bonding  model  2 

The  above  discussions  are  based  on  the  fact  that  only  the  inclusion  is  subject  to  an  axial 
load  while  the  plate  is  free  of  any  external  loads.  Such  a  model  is  often  used  to  illustrate 
how  the  load  is  being  transferred  from  the  fiber  to  the  matrix,  as  we  have  reviewed  in  the 
Introduction.  But  in  reality  in  composites,  both  the  fiber  and  the  matrix  carry  portions  of 
the  external  load  provided  that  the  shear  modulus  p.i  of  the  matrix  is  not  negligible.  For  this 
reason  in  Fig.  3  we  consider  a  modified  model  which  allows  the  plate  to  carry  an  axial  load 
Gi.  The  magnitude  of  g,  may  be  different  from  that  of  Gq  and  in  our  subsequent  numerical 
calculations,  the  ratio  Gj  /Gq  will  be  considered. 
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Most  of  the  formulas  developed  for  perfect  bonding  model  1  still  hold  except  some 
modiHcations  need  to  be  made  for  expressions  (2.15-2.17). 


u<pXi)^. - 2 - r 

»  2pi(mi+l) 

(2.31) 

(2.32) 

„(rta)  _55ia_. 

(2.33) 

oS:x'>^»'=0 

(2.34) 

OS’"’-®! 

(2.35) 

^XD^XD^Xl)^ 

(2.36) 

On  the  other  hand,  the  particular  solution  for  the  inclusion  is  the  same  as  that  given  in 
(2.18-2.23). 

Imperfect  Bonding  Model 

Perfect  bonding  conditions  require  that  the  displacements  and  the  stresses  be  continuous  at 
the  interface  throughout  the  thickness.  Since  both  of  the  above  models  satisfy  (2.9)  to 
(2.10),  we  define  them  as  perfect  bonding  cases.  As  it  was  previously  noted,  imperfect 
bonding  models  which  have  been  examined  based  on  2D  consideration  are  basically  of  two 
types:  in  the  first  type  the  transfer  of  load  across  the  interface  is  described  by  (Coulomb 
friction,  while  in  the  second  the  transfer  of  load  across  the  interface  is  described  by  a 
continuous  elastic  spring  connecting  the  fiber  and  the  matrix.  In  this  analysis  we  choose  the 
latter  for  it  is  more  suitable  to  the  structure  of  our  complementary  solution  given  in  the 
subsequent  analysis. 

This  model  is  shown  schematically  in  Fig.  4.  Oucial  to  the  problem  of  interest  is  the 
bonding  between  the  fiber  and  the  matrix  at  r=a  where  the  relevant  stresses  as  well  as  the 
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displacements  u  and  v  must  be  continuous.  The  vertical  displacement  w,  however,  is 
allowed  to  be  discontinuous  but  its  jump  must  be  proportional  to  the  shear  stress,  i.e.. 


along  p=a: 


u">-u®=K.  £>;  (i.lor2) 


where  1/Ks  is  the  interface  shear  stiffness. 


(2.37) 


A  boundary  condition  of  the  type  (2.37)  is  inappropriate  for  the  edge  region  (near  z=h)  for 
the  displacement  field  is  non-singular  while  the  stress  field  is  singular  (Folias,  1989). 
However,  except  in  the  boundary  layer  region,  the  results  arc  expected  to  be  valid. 

In  view  of  this  consideration,  the  last  boundary  condition  in  (2.14)  becomes 


u»'.u®=K.t®:  (i=lor2) 


(2.38) 


All  other  boundary  conditions  remain  the  same.  Finally  in  terms  of  the  complementary  and 
particular  solutions,  the  boundary  condition  (2.38)  becomes 


All  other  boundary  conditions  are  given  by  eq.  (2.24-2.30). 


(2.39) 


It  remains,  therefore,  for  us  to  find  a  complementary  solution  such  that  it  satisfies  Navier's 
equations  and  the  appropriate  boundary  conditions. 


The  general  method  for  constructing  solutions  to  some  3D  mixed  boundary-value  problems 
has  been  developed  by  Folias  (1975, 1990).  The  solution  may  be  written  in  the  following 
more  convenient  form  by  considering  6-independence  of  the  problem: 


(i) 


I 

u(®W>=^^^I^-^^{2(mi-l)fi(Pvz)+niif20vz)) 

j  82X5> 

^  y  dx  ^mi+1  ^^dxSy 

oo 

^  ^  1  (Pv2)+mif20vz)  1 

i  (i)  +  1  (i).  V  ^  —L.  ,2_. 

^mi+1  ^3  2  y  mi+1  *  9x2 


(3.1) 


an 


(i) 


fO 


Jj-jf  l<"'i-2)f3(M-nii£4(3vz))-  ^  ^ 

Furthermore  the  stresses  are  given  by  the  Hooke's  law; 


(3.2) 


(3.3) 


(i) 


^X{2Pv-^fl(Pvz)  +  ~^2(mi-l)fi(Pv2)  +inif2(Pvz)]) 

(3.4) 


fO 


dxf  a2x§^  2  1 

'  ^x  ‘^ax^  *  mi+l  dy  'mi+l^^ax^ay 


J{2pJ-5r-fl(M)  -( "5^  -  ^-^)P(tni-l)fl(Pvz)  +mife(Pvz)]) 


(3.5) 


<  2ini  <  1 

"  dx  ^  9x^  ^  mi+l  "SbT '  mi+l  ^^9x29y 


3^’  ^  n^-l  3^?  I  ,3’^" 

*  dx  ^  ^dxSy  '*’  mj+l  3x  '  mj+l  dx? 


(3.7) 


1  m;  •*  ^  ^  t/ 

^  ^  ^ii-S57^v{f3(M)  +  f4(Pvz)) 

1  m;  ~  ^^^2^ 

2h<’'‘’=  ■  J,-^ff3(l>v*)  +  «4(M)) 


where 


ctn=  L  t  n— l,2,3f.... 


Py  are  the  roots  of  the  equation 


(3.10) 


sin(2Pvlt)s  ~  2pvlt 


(3.11) 


(i=l,  2)  are  functions  of  x  and  y  which  satisfy  the  reduced  wave  equation: 


(^^■I«)hS’-o 


(3.12) 


kf\X^^>andX§^  are  two  dimensional  harmonic  functions,  and 


f]  Ovz)’’cos0vh)cos(^z) 

f2(Pvz)»^hsin(Pvh)cos(Pvz)  -  ^zcos(Pvh)sinOvz) 


(3.13) 

(3.14) 


f3(Pvz)^osOvh)sinOvz) 

f4(^z)=^hsin(pvh)sin(Pvz)  +  ^vzcos(^h)cos(^z) 


(3.15) 

(3.16) 
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The  reader  may  notice  that  the  above  coixq)iementaiy  solution  automatically  satisfies  the 
boundary  conditions  at  Izl^h  and  that  it  only  remains  to  satisfy  the  remaining  boundary 
conditions  at  the  interface  r=a. 

In  the  case  of  perfect  btmding  (model  1  and  model  2),  we  transform  the  complemmtary 
solution  into  cylindrical  coordinates  in  order  m  talm  advantage  of  the  symmetry.  Thus, 


equations  (2.25*2.30)  may  be  written  in  the  following  form: 
su.2e«4'">  cos^e  «in(2e)(x<«’> 

jsiitzeco'J"’  ^4’®V|siii2e  (o'Jjf' ’  4cosaexx‘^‘’  =o  p.  1 8) 

sin0(t^;^*^  -x^^=0  (3.19) 

sin0(u(«Xi).u(cX2))  +  cose(v('Xi)-v(cX2))=Ci a  +  z  -7^  a  (3.20) 

2pi(mi+l) 

cos0(u<®Xl)-u<®X2)) .  sin0(v(cXi).v(cX2))s^  (3.21) 


(3.22) 
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In  as  much  as  model  1  is  a  special  case  of  model  2  wth  ai=0.  from  here  on  we  will 
concentrate  on  the  recovery  of  the  solution  to  model  2. 


Taking  into  account  the  symmetry  of  the  problem  (6  independence),  we  let 


^1) 

a  CivKoOvt)  (3.24) 

*  C2vIo(^^*‘)  (3.25) 

and  all  X«0  excqrt 

X?^=ysine  (3.26) 

(3.27) 


where  Iq  and  Ko  are,  respectively,  the  modified  Bessel  functions  of  the  first  and  second 
land  of  order  zero;  Civ  and  C2v  ut  ctmiplex  constants.  The  constants  A  and  Ci  are  to  be 
determined  in  such  a  way  that  all  remaining  conditions  are  satisfied. 


In  the  case  of  imperfect  bonding,  the  boundary  condition  corresponding  to  (3.22)  becomes 


w(cX1).w(‘=X2)=k,x®-».[ 

a 


Op  1-2v2  2v2  p  ■,  miOi 
2P2  1*V2  ’1.V2  "^■2pi(mi+l)* 


(3.28) 


while  all  other  boundary  conditions  are  obtained  from  (3.17-3.21)  by  letting  Oi  =0. 


Perfect  Bonding  Model  1 

In  as  much  as  the  interface  is  of  greatest  practical  interest,  plots  at  r=a  are  given  for  the 
displacements  and  stresses  as  functions  of  z.  The  case  of  perfect  bonding  (model  1)  will  be 
taken  as  the  basis  of  our  first  discussion.  Two  different  shear  moduli  ratios  are  chosen: 
these  are  and  li2/lti=17  corresponding  to  (a)  and  (b)  respectively  in  each  Figure. 

The  purpose  of  these  plots  is  to  examine  how  the  shear  moduli  ratio  affects  the 
displacement  and  the  stress  fields.  Each  plot  consists  of  several  curves  which  corresponds 
to  different  a/h  ratios.  Superscripts  are  used  to  distinguish  whether  '’he  stress  or 
displacement  is  in  the  fiber  or  in  the  matrix  as  they  have  been  defined  early  in  the  Figure  1. 
However  the  superscripts  are  omitted  if  the  stress  or  displacement  is  same  for  the  fiber  and 
the  matrix  at  the  interface. 

Stress  Off  (Fig.5  (a)  (b) )  increases  with  z/h.  At  the  center  region,  the  Figure  shows  only  a 
slight  change  of  the  magnitude  of  the  stress.  However,  it  jumps  suddenly  to  a  large  value 
near  the  plate  surface,  suggesting,  therefore,  the  presence  of  a  singularity  (Folias,  1989).  It 
is  interesting  to  note  that  in  the  cases  of  small  a/h  ratios,  is  slightly  negative  over  a  wide 
span  of  the  center  region.  This  suggests,  therefore,  that  small  compression  exists  at  the 
interface  between  the  fiber  and  the  matrix  at  this  region.  In  the  cases  of  large  a/h  ratios, 
however,  the  positive  value  of  a„  suggests  that  small  tension  exists  at  the  interface.  Our 
analysis  shows  that  the  shear  stress  Xn  (see  Fig.  6  (a)  (b))  starts  from  zero  at  z^O  and  then 
increases  in  absolute  value  along  the  interface.  A  sudden  jump  of  the  stress  near  the  points 
zsh  is  observed  and  may  very  well  explain  why  a  crack  is  most  likely  to  initiate  at  this 
region.  The  sign  for  the  shear  stress  is  negative  because  the  shear  stress  will  always  be 
opposing  the  direction  of  the  external  stress  Cq  applied  at  the  face  of  die  fiber.  The  variation 
of  the  magnitude  of  tire  shear  stress  Xn  seems  oaore  complicated  in  the  cases  with  p2/l^i=17 
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than  in  the  cases  with  because  the  fonner  one  (Fig.  6  (b))  shows  the  shear  stress 

to  lecum  as  aAt  increases  after  it  reaches  a  maximum  absolute  value  at  a  certain  a/h  value.  It 

t^pears  that  both  *Very  thick"  and  “very  thin"  plates  will  induce  very  small  shear  stress  at 

(2) 

the  center  region.  The  matrix  tensile  stress  at  the  interface  is  plotted  in  Fig.  7(a)  (b), 

the  stress  increases  as  ratio  a/h  increases.  A  extremely  high  stress  in  the  vicinity  of  fiber 

may  cause  the  matrix  to  break  at  this  region.  For  most  cases,  the  stress  shown  in 
p^g  8  (a)  (b)  decays  as  z/h  moves  from  the  fiber  end  (z=h)  to  the  fiber  center  (z^O).  But 

for  the  cases  with  a/h=1.0,  the  result  shows  that  the  stress  bounces  back  slightly  when  it 
reaches  the  fiber  center.  It  is  also  noted  diat  the  tensile  stress  at  the  matrix  will  decrease  as 
increases  by  comparing  Fig.  7(b)  with  Fig.  7(a)  while  the  tensile  stress  at  the  fiber 
will  increase  by  comparing  Fig.  8(b)  with  Fig.  8(a).  This  suggests  that  the  load  diffusion 
ftom  the  fiber  to  the  matrix  will  decrease  when  the  fiber  is  stiffer. 

The  displacement  u^  at  z=h  vs  the  radial  direction  is  given  in  Fig.  9.  It  shows  that  a  softer 
fiber  (pi/^i^^)  gives  a  smoother  conn^on  of  the  displacement  at  r=a  while  a  stiffer  fiber 
(P2Ati=17)  gives  a  sharper  connection  of  die  displacement  at  the  interface. 

Perfect  Bonding  Model  2 

In  as  much  as,  in  practice,  the  matrix  also  carries  a  small  portion  of  the  applied  load  it  is 
desirable  now  to  extend  our  model  to  the  case  where  die  matrix  also  carries  part  of  the  load 
in  the  direction  of  the  fiber  length.  Thus,  perfect  bonding  noodel  2  allows  us  to  examine 
how  the  load  is  being  transferred  ftom  the  fiber  to  the  matrix  and  vice  versa.  As  it  was 
previously  noted,  the  perfect  bonding  model  1  is  taken  as  the  basis  of  our  discussion.  To 
obtain  the  effect  of  Oj  (recall  Fig.3),  the  shear  moduli  ratio  is  chosen  as  p2/lti=17  and  the 
geometric  parameter  a/h  is  fixed  at  a/h=0.05.  Three  curves  are  plotted  in  each  figure,  which 
corresponds  to  Oi/ao=0.0,  0.05  and  0.1,  respectively.  It  is  worth  pointing  out  that  ratio 


Oi/ao  is  not  an  independent  parameter.  Actually  this  ratio  will  depend  on  the  shear  moduli 
ratio  |i2A^i  if  a  uniform  strain  is  assumed  on  the  surface  of  fiber  and  matrix.  Thus,  for  the 
parameters  chosen,  a  realistic  ratio  of  o/obmust  be  around  1/17. 

Stress  Git  shown  in  Fig.  10  will  not  be  greatly  influenced  in  the  center  region  as  Oj/Co 

increases.  However  it  may  jump  to  the  negative  direction  when  Oi/Oo  passes  a  certain 

critical  value.  The  numerical  calculations  appears  suggesting  that  the  stresses  (including 

other  stresses  which  will  be  discussed  later)  will  change  the  direction  of  the  singularity 

when  Ci/ao>pi/42.  Such  a  change  indicates  that  the  debonding  in  the  edge  region  becomes 

less  possible  because  of  the  negative  value  of  the  radial  stress.  Fig.  1 1  shows  the  shear 

stress  trz  vs  z/h.  Identical  to  an  increase  of  gJcq  beyond  a  certain  value  will  also 

change  the  sign  of  t^,  which  implies  that  the  shear  stress  now  will  be  in  the  same  direction 

of  external  tensile  stress.  Moreover,  this  figure  suggests  the  presence  of  a  negligible  shear 

C2) 

stress  at  the  interface  when  Oj/ao,  ii2/^iand  a/h  are  properly  chosen.  Fig.  12  shows  o  jz 
the  interface.  Notice  that  the  magnitude  of  the  stress  dramatically  increases  as  the  ratio  of 
oi/oo  increases  (  see  the  curve  with  Gi/GO=0.1).  This  figure  also  indicates  that  the  stress 
value  at  die  center  (z^O)  may  be  the  maximum  value. 

Imperfect  Bonding  Case 

When  the  interface  is  allowed  to  slip,  the  magnitude  of  the  stresses  in  the  matrix  will  be 
expected  to  decrease  while  in  the  fiber  they  are  expected  to  increase  slightly.  In  order  to 
show  the  effect  of  stiffness  of  the  interface,  two  parameters  are  fixed:  they  are  a/h=0.05 
and  H2/p.i=l7.  Three  curves,  which  corresponds  to  the  reciprocal  of  the  spring  constants 
Ke=0,0, 0.1  and  1.0,  are  plotted. 


20 


The  radial  stress  o„  shown  in  Fig.  13  increases  slightly  as  the  interface  becomes  softer 

(note  Kes0,0  corresponds  to  perfect  bonding  case,  or,  the  rigid  interface).  This  indicates 

that  the  softer  interface  will  increase  the  chance  of  inteifacial  debonding.  The  shear  stress 

Xrz  shown  in  Hg.  14,  on  the  other  hand,  shows  a  dramatic  decrease  of  the  stress  along  the 

interface.  It  appears  that  the  interface  is  somehow  relaxed  by  a  softer  interface.  The  fiber 
(2) 

tensile  stress  (Fig.  15)  increases  as  slippage  is  allowed  to  increase.  This  suggests  that 
less  load  will  be  transferred  firom  the  fiber  to  the  matrix  if  the  interface  is  softer,  which  also 
substantiates  our  initial  assumption. 

The  spring  constant  can  be  calculated  by  the  equations  given  by  Stief  and  Hoysan  (1986), 
where  the  coated  fiber  model  and  the  interfacial  crack  model  are  considered.  As  pointed  out 
by  them,  a  significant  change  of  stress  Held  can  only  happen  when  the  coating  material  is 
very  soft.  This  means  that  the  stress  field  calculated  from  the  practical  parameters  may  be 
quite  close  to  that  calculated  from  the  perfea  bonding  model. 
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COMPARISONS  AND  CONCLUSIONS 

Our  present  results  are  consistent  with  the  previous  observations  based  on  ID  and  2D 

(2) 

considerations  (Lawrence,  1972  and  Banbaji,  1988)  that  Xn  and  will  attain  their 

maximum  values  at  the  loaded  end  of  the  fiber.  Although,  our  present  model  differs  fiom 

previous  3D  models  ( Muki  et  al.  1970,  Luk  et  al.  1979,  Haener  et  al.  1967),  the  stress 

profiles  obtained  are  similar.  For  exan:q)le,  conqiarin  g  die  result  of  Haener  (1965)  to  that  of 

(i)  (i) 

our  perfect  bonding  model  1,  one  finds  that  the  stresses  and  o,  obtained  in  these  two 

different  models  both  show  a  flat  behavior  in  the  center  region  (acmally  the  value  is  rather 

small  in  this  region  for  a„ )  and  then  a  sudden  jump  to  a  larger  value  near  the  surface.  In 
addition,  our  present  numerical  results  do  confirm  the  presence  of  a  sness  singularity. 

The  qualitative  explanation  of  fiber  multiple  crack  phenomenon  found  in  the  Dogbone  test 
samples  (DiBenedetto  et  al.,  1986  and  Bascom  et  al.,  1986)  can  be  made  by  investigating 
the  fiber  tensile  stress  along  its  length.  Fig.  16  shows  two  curves  corresponding  to  two 
different  geometric  ratios  a/h=0.01  and  a/haO.OS,  where  [a]  denotes  the  tensile  strength  of 
the  fiber.  As  we  can  see,  the  tensile  strength  in  longer  fiber  (a^=0.01)  has  passed  the  dash 
line  and  hence  will  break  while  the  shorter  fiber  (a/h=O.OS)  will  not  break  because  the 
maximum  tensile  stress  is  below  the  dash  line.  However,  if  Oj  increases,  the  dash  line  will 
move  lower  provided  that  the  tensile  strength  of  the  fiber  is  independent  of  the  ratio  a/h. 
This  may  lead  to  the  breakage  of  shorter  fiber.  It  becomes  clear  that  the  fiber  tends  to  relax 
itself  by  having  shorter  length.  Another  interesting  phenomenon  is  that  the  non- 
dimensionalized  tensile  stress  at  the  fiber  center  for  the  case  a/h=0.01  is  about  17,  which  is 
equal  to  the  shear  moduli  ratio  of  the  current  problem.  Cox  (1952)  has  predicted  that  the 
maximum  tensile  stress  of  the  fiber  under  this  model  will  occur  at  the  fiber  center  and  the 
value  wiU  be  E20j/Ei.  This  will  yield  170i  assuming  the  Poisson  ratio  is  same  for  the  fiber 
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and  the  matrix  in  our  current  problem.  The  prediction,  however,  becomes  inaccurate  when 
the  fiber  length  decreases.  This  indicates  that  the  E20i/Ei  can  only  be  taken  as  a  uplimit  of 
the  fiber  tensile  stress. 

A  quantitative  comparison  of  our  current  results  with  the  experimental  results  obtained  by 
Tyson  and  Davies  (1965)  is  given  in  Fig.  17.  The  parameters  of  their  2D  experimental 
model  are  taken  to  fit  our  current  3D  model.  The  diameter  of  the  fiber  is  tal»n  as  4nim  and 
the  half  thickness  of  the  plate  h  is  determined  by  photoelastic  method  by  measuring  the 
distance  from  the  fiber  end  to  the  isotropic  point  All  other  parameters  are  given  in  the  paper 
and  can  be  used  directly.  As  seen  from  the  Figure,  our  result  predicts  the  interfacial  shear 
stress  quite  satisfactorily  throughout  the  interface.  A  deviation  starts  at  x=4mm,  one  fiber 
diameter  from  the  fiber  end.  The  authors  are  aware  that  this  deviation  may  be  caused  by 
neglecting  the  stress  concentration  effect  at  the  edge  region  as  discussed  previously.  On  the 
other  hand,  it  also  indicates  that  the  influence  of  such  a  stress  concentration  effect  is 
extremely  localized.  In  addition,  the  result  calculated  from  the  shear  lag  analysis  greatly 
underestimates  the  interfacial  shear  stress  especially  in  the  vicinity  of  the  fiber. 

In  view  of  the  numerical  results,  one  may  conclude  that: 

(1)  The  thickness  as  well  as  the  material  properties,  play  a  fundamental  role  on  the  failure 
mechanism  of  a  plate  with  a  cylindrical  inclusion.  More  specifically,  the  geometric 
parameter  a/h  and  the  shear  moduli  ratio  are  two  important  ratios  which  greatly  affect 
the  displacement  and  stress  fields. 

(2)  A  boundary  layer  effect  is  shown  to  exist  in  the  vicinity  of  z=h  where  the  presence  of  a 
stress  singularity  (Folias,  1989)  may  result  to  crack  initiation. 

(3)  The  shear  lag  analysis  may  underestimate  the  interfacial  shear  stress  and  may 
overestimate  the  fiber  tensile  stress. 
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(4)  When  Oi/ao<kLi/M-2>  the  interfacial  debonding,  slippage  and  fiber  breakage  will  most 
likely  initiate  at  the  edge  region. 

(5)  When  Gi/ao>Hi/^i2,  the  interfacial  slippage  will  initiate  at  the  edge  region  while  the 
interfacial  debonding  and  fiber  breakage  may  initiate  at  the  center  region. 

(6)  Taking  a  uniform  tensile  stress  at  the  the  plate  surface  by  neglecting  the  stress 
concentration  effect  at  the  edge  re^on  may  slightly  underestimate  the  interfacial  shear  stress 
near  the  fiber  end.  However,  such  a  error  reduces  quickly  as  the  distance  from  the  fiber  end 
increases. 

(7)  The  substitution  of  the  "interphase"  with  an  elastic  spring  may  give  inaccurate 
information  of  the  stress  at  the  edge  region  because  one  of  the  boundary  conditions  breaks 
down  at  the  re^on. 
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Fig.  1  Perfect  bonding  Model  1:  only  fiber  is  subjected  to  an  axial  load 


Fig.  3  Perfect  bonding  Model  2;  both  fiber  and  matrix  are  subjected  to  axial  loads 


(b)  H24ii=17  and  vi=v2=0.33 
Fig,7  Stress  vs  z/h  for  perfect  bonding  Model  1 
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16  Stress  O^Vs  z/h  for  short  fiber  model  where  VjaV2»0.33;  Pj  f  ^7 


Fig.  17  G)niparison  with  the  experimental  results  abtained  by  Tyson  and  Davies 
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The  previous  analysis  may  now  be  extended  to  the  the  case  of  a  periodic  array  of 
fibers  embedded  into  a  matrix  plate  and  under  the  action  of  a  uniform  load  in  a  direction 
parralel  to  the  fiber  axis  (see  figs  1  and  2) .  Without  going  into  the  mathematical  details, 
which  are  similar  to  the  previous  parts,  the  results  are  given  in  figures  3  through  11.  From 
these  rigures  it  becomes  evident  that  the  load  transfer  characteristics  from  the  fibers  to  the 
matrix  are  gavemed  by  the  material  properties,  the  fiber  radious  to  fiber  lentgh  ratio,  the 
fiber  volume  fraction,  and  finally  the  load  ratio  matrix  and  fibers. 


1.  Geometrical  configuration 


Fig.  SL  Perfect  bonding  Model  both  fiber  and  matrix  are  subjected  to  axial  loads 


sitas45  z<^0 


G2/G1  =  16.67  v1s0.34  v2«0.22 


Fig.  4  The  stress  Ozz  at  the  interface  versus  a/b  at  2=0 , 0=0  and 
fl2=  Oi/oo=0.1. 
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Fig.  10  The  stress  at  the  interface  versus  z=h  at  8=  45  and  fl2: 
Oi/oo=0.1  and  various  a/h  ratios. 
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ABSTRACT 


The  adhesive  butt  joint  has  been  so’v^ed  by  using  the  theory  of 
linear  elasticity  under  the  assumption  of  plane  strain  conditions.  Both 
the  adherend  and  the  adhesive  layer  are  assumed  to  be  homogeneous 
and  isotropic  materials  with  different  material  properties.  Perfect 
bonding  is  also  assumed  to  prevail  at  the  interface. 

The  numerical  results  show  that  under  a  tension  load  the 
interface  plane  does  not  remain  plane  and  that  the  normal  stress 
possesses  a  weak  singularity  near  the  free  surface  z/h=l.  It  is  also 
found  that  there  is  significant  variation  in  the  normal  and  shear  stress 
profiles  across  the  adhesive  thickness  especially  near  the  free  surface 
and  that  the  stress  distribution  is  very  sensitive  to  the  value  of  the 
adhesive  thickness  to  width  ratio  (a/h)  particularly  when  this  is  less 
than  1.  An  examination  of  the  octahedr£d  shear  stress  shows  that  a 
crack  is  most  likely  to  initiate  at  the  comer  of  the  adhesive  layer. 
Finally,  the  numerical  results  show  a  good  agreement  with  the  results 
obtained  by  finite  element  method  (FEM). 


INTRODUCTION 

In  recent  years,  there  has  been  a  continuing  increase  in  the  use 
of  adhesives  as  an  alternate  method  for  assembling  the  structures. 
From  the  definitions  of  ASTM  (1982),  an  adhesive  is  a  substance 
capable  of  holding  materials  together  by  surface  attachment,  an 


adherend  is  a  body  that  is  held  to  another  body  by  an  adhesive  cind 
adhesive  joint  is  the  location  at  which  two  adherends  are  held 
together  with  a  layer  of  adhesive.  Figure  1  shows  some  adhesive  Joints 
that  are  commonly  found  in  engineering  practice. 

The  single-lap  joint  is  one  of  the  most  commonly  found  joint 
configurations  in  engineering  practice  (see  Figure  1).  There  are 
numerous  papers  in  the  literature  that  deal  with  the  problem  of  a  lap 
joint  from  a  theoretical  as  well  as  a  finite  element  point  of  view. 
However,  very  few  analyses  have  been  done  on  adhesive  butt  joints. 
Adhesive  butt  joints  are  usually  designed  with  rectangular  or  circular 
cross-section.  The  mechanical  properties  of  adhesive  in  bulk  form  are 
different  from  those  in  thin-film  form,  so  the  adhesive  butt  joint  has 
been  widely  used  as  a  testing  specimen  to  obtain  the  thin-film 
mechanical  properties  of  adhesive.  Gent  et  al.  (1970)  have  analyzed 
the  incompressible  adhesive  (i.e.,  its  Poison’s  ratio  is  0.5)  in  tension  or 
compression  between  rigid  adherends.  Lindsey  (1966)  has  calculated 
the  stress  distribution  of  the  adhesive  butt  joint  in  tension  with  the 
assumptions  that  adherends  are  rigid  and  stresses  are  constant  across 
the  adhesive  layer.  Two-dimensional  plane  stress  finite  element 
analysis  has  been  done  by  Harrison  et  al.  (1972)  who  use  the  same 
assumptions  as  Lindsey  (1966)  and  maximum  shear  stress  is  found 
near  the  free  surface. ,Alw£ir  et  al.  (1976)  and  Adams  et  al.  (1978)  have 
used  the  finite  element  method  to  analyze  plane  stress  and 
axisymmetric  adhesive  butt  joint.  Both  anadyses  assume  that  the 
adherends  are  nonrigid  materials  and  their  numerical  results  show 


that  there  is  significant  variation  in  normal  and  shear  stress  profiles 
across  the  adhesive  thickness.  A  similar  conclusion  is  obtained  by  Sawa 
et  al.  (1989)  who  analyze  the  axisymmetric  adhesive  butt  joint  under 
tension  by  using  the  Michell's  stress  function  from  the  three- 
dimensional  theory  of  elasticity. 

In  this  paper,  a  plane  strain  adhesive  butt  joint  loaded  in  tension 
has  been  anal5rzed  using  the  general  solution  derived  by  Folias  (1975, 
1990).  Both  the  adherend  and  the  adhesive  layer  are  assumed  to  be 
homogeneous  and  isotropic  materials  with  different  material 
properties.  Perfect  bonding  (i.e.,  continuity  of  displacements  at  the 
interface  between  adhesive  and  adherend)  is  assumed  to  prevail  at  the 
interface.  The  stress  distribution  of  adhesive  layer  with  respect  to 
different  values  of  geometric  ratio  (a/h;  adhesive  thickness/  adhesive 
width)  and  material  property  ratio  (G2/G1;  adhesive  shear  modulus/ 
adherend  shear  modulus)  have  been  investigated  and  these  results  give 
further  understaindlng  of  the  adhesive  butt  joint  problem  in  the 
presence  of  a  tensile  load. 


FORMULATION  OP  THE  PROBLEM 

Consider  two  adherends  of  thickness  2h  bonded  together  in  a 
butt  joint  configuration  by  an  adhesive  layer,  as  shown  in  Figure  2.  The 
adherends  occupy  the  space  a  ^  Ixl  <  «>,  lyl  <  «>,  Izl  <  h  and  the 
adhesive  layer  occupies  the  space  lx  I  ^  a,  lyl  <  «>  and  Izl  <  h.  Let  the 
joint  be  subjected  to  a  uniform  tensile  load  a©  along  the  x-axis  and 


parallel  to  the  bonding  planes.  The  shear  modulus  G2  and  Poisson’s 
ratio  V2  represent  the  material  constants  of  the  adhesive  layer.  In 
order  to  simplify  this  adhesive  butt  joint  into  a  symmetrical  problem, 
one  may  assume  that  these  two  adherends  are  homogeneous,  isotropic 
and  linearly  elastic  materials  with  the  same  shear  modulus.  Gj,  and 
Poisson's  ratio,  v^,  respectively.  Since  this  adhesi\«j  butt  joint  has  been 

simplified  to  a  symmetrical  problem,  one  needs  only  to  consider  the 
one  quarter  part  of  this  joint  with  displacement  u  equal  to  zero  along 
the  line  x  =  0  and  displacement  w  equal  to  zero  along  the  line  z  =  0. 
According  to  Figure  2.  this  adhesive  butt  joint  has  a  very  large 
dimension  in  the  y  direction  and  it  is  reasonable  to  analyze  it  as  a  two- 
dimensional  plane  sir?in  problem. 

For  this  type  of  problem,  it  is  convenient  to  separate  the  solution 
in  two  parts,  the  particular  and  complementary  solution 


u 


ID 
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1  ) 


=  u‘P’“’  +  u'^»‘' 

=  w'  P  +  w‘  "  »  ‘  ’ 


(1) 

(2) 


i=1.2.  where  i=l  represents  the  adherend  and  i=2  represents  the 
adhesive  layer.  Superscript  (p)  represents  the  particular  solution  and 
superscript  (c)  represents  the  complementary  solution. 


The  particular  part  must  satisfy  the  Navler’s  equations  and  the 
loading  conditions  which  without  the  presence  of  the  adhesive 


bonding  interface,  and  the  complementary  part  takes  care  of  the 
adhesive  bonding  interface  with  no  contribution  far  away  from  the 
interface.  The  particular  solution  of  this  problem  can  be  constructed 
as: 


(p)(l)  "'i'l 
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2G. 

(3) 
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(4) 

lplli)_ 
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“  -  m, 
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(5) 

-S'"  =  -0 

(6) 
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(7) 

-S'"  =  0 

(8) 

Assuming  that  a  particular  solution  has  been  constructed,  one 
then  needs  to  find  four  complementary  displacement  functions 
u(c)(i)(x,z),  w*‘^**^*(x,z),  u^^**^*(x,z)  and  w^^^^^^^lx.z)  in  order  to  satisfy  the 
Navier's  equations  and  the  required  boundary  conditions  which  are: 

ixi^oo,  =  0 

Izl  =h,  0 


(10) 


v/r 


(0(2)  _  (0(2)  _ 

^zz  ^xz 

(11) 

(p)(l)  .  _(C)(1)  _  _(p)(2)  (c)(2) 

■*'  -  ^XX  +  ^XX 

(12) 
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^xz  ^xz  '^xz  ^xz 

(13) 

^(p)(l)  ^  ^(0(1)  ^  ^(P)(2)  ^  ^(0(2) 

(14) 

y(p)(l)  ^  ^(0(1)  _  y(pK2)  ^  ^(0(2) 

(15) 

Since  the  adherend  is  a  semi-infinite  strip,  the  complementary 
displacements  and  stresses  of  the  adherend  must  vanish  as  x  tends  to 
infinity.  On  the  other  hand,  the  complementary  displacements  and 
stresses  of  the  adhesive  layer  must  also  be  bounded  at  x=a. 

METHOD  OP  SOLUTION 

A  general  three-dimensional  solution  to  Navier's  equation  for 
plates  of  uniform  thickness,  2h,  and  with  plate  faces  free  of  stress  has 
been  constructed  by  Folias  (1975)  eind  applied  to  the  problem  of  a 
plate  of  finite  thickness  containing  a  finite,  throu^  the  thickness,  line 
crack.  Later,  Folias  et  al.  (1990)  applied  the  solution  to  the  problem  of 
a  finite  thickness  plate  that  has  been  weakened  by  a  cylindrical  hole. 
Without  going  into  the  mathematical  details,  the  general  form  of 
solution  obtained  by  Folias  et  al.  (1990)  is: 
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Vsl 
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Since  the  adhesive  butt  joint  considered  here  is  a  two- 
dimensional  plane  strain  problem,  one  may  simplify  the 
complementary  displacements  as: 
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(20) 


]£  \  {  I  (m.-2)  cos(p^h)  -  m,p^h  sin(p^h)  ] 

“  1  ^  V=1 

sin(p^z)  -  nijp^z  cos(p^h)  cosCP^z)  } 


By  using  the  Hooke’s  law.  one  may  find  the  corresponding 
complementary  stresses  as: 
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where  i  =  1,  2,  Ij^^*  =  0,  and  Py  are  the  roots  of  the  transcendent2d 
equation 


sin  (2pyh)  =  -  (2Pyh), 


(24) 


It  is  easy  to  assure  that  the  complementary  stresses,  i.e.,  Txz**^’*** 
,vnd  a^z*^^****.  vanish  at  the  plate  faces  lzl=h. 

Substituting  equations  (19)  and  (20)  into  the  Navier's  equations  , 
one  can  show  that  the  complementary  displacement  field  does  indeed 
satisfy  Navier's  equations,  provided  that  the  function  Hy***  satisfies  the 

following  equations: 


(25) 


The  proper  type  of  functions  and  for  this  problem  is: 


.(2) 


(x)  =  cosh(  PyX  ) 


(26) 

(27) 


By  substitution  the  particular  and  complementary  solutions  into 
the  guided  boumdary  conditions  (12)-(15),  one  may  conclude  that  the 
governing  boundary  conditions  become  respectively: 
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where  and  A^*^*  are  the  normalized  coefficients,  with  the 

definitions  of: 
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Examination  of  equations  (28)-(33)  shows  that  the  coefficients 
Ay*^*  and  A^*^’  are  functions  of  Py  ,  the  dimensionless  ratio  a/h,  the 

shear  modulus  ratio  Gj/Gi  and  the  Poissons  ratio  V1/V2.  Since  the 

roots  p2’  P4-  P6 .  sure  the  complex  conjugates  of  p^.  P3.  ps .  one 

concludes  that  the  unknown  complex  coefficients  A2**’.  Ai'*’,  Ag*** . 

are  also  the  complex  conjugates  of  Aj'**,  Ag***,  Ag**' .  (where  i  =  1,2). 

The  unknown  complex  coefficients  Ay^**  and  Ay*^*  are  to  be 
determined  from  the  system  of  linear  equations  (28)-(31),  by  using 
the  method  discussed  in  Kantorovich  and  Krylov  (1964,  p. 54-56).  The 
system,  equations  (28)-(31).  is  veiy  sensitive  to  even  small  changes  of 
the  coefficients,  so  the  methods  of  "collocation"  and  "least  squares" 
lead  to  a  nonconvergent  solution.  However,  the  method,  which  was 
discussed  by  Kantorovich  and  Krylov  (1964),  shows  the  complex 
coefficients  to  converge  as  the  number  of  roots  taken  increases. 


Once  the  unknown  complex  coefficients  have  been  determined, 
the  displacement  and  stress  fields  may  then  be  computed  at  any  point 
in  the  adhesive  layer. 

The  displacement  field  in  the  adhesive  layer  becomes: 

(  ^)  U-  =  iVi  X  ,  £  cos(p„h, 

2m^  ^  cosh(p^a)  ^  2 
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Similarly,  the  stress  field  of  the  adhesive  layer  becomes: 
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NUMERICAL  RESULTS  AND  DISCUSSION 

In  the  present  analysis,  we  let  the  a/h  =  0.1.  G2/G1  =  0.05.  = 

0.3  and  V2  =  0.33.  After  determining  the  complex  coefficients  and 
Av*^’,  one  needs  to  substitute  these  complex  coefficients  into  equations 
(28)  -  (31)  to  examine  how  well  the  complex  coefficients  satisfy  the 
original  boundary  equations  (28)  -  (31).  Figure  3  shows  how  the  error 
of  equation  (28)  improves  as  one  increases  the  number  of  roots  from 
100  to  200.  The  numerical  study  shows  that  the  errors  of  these 
boundary  equations  improve  only  3%  if  one  increases  the  root  number 
V  from  200  to  220  and  it  is  reasonable  to  set  200  roots  throughout  this 

present  euialysis.  If  further  accuracy  is  required,  then  more 
sophisticated  algorithms  must  be  used  for  the  solution  of  the  matrix 
system. 

By  substituting  the  complex  coefficients  Ay*^*  into  equations  (34) 
-  (38).  one  obtains  the  displacement  and  stress  fields  in  the  adhesive 
layer.  The  numerical  study  indicates  the  interface  does  not  remain 
plane  under  load,  although  the  shear  modulus  of  adherends  is  20 
times  as  stiff  as  that  of  adhesive.  The  variation  of  displacement  u*^' 
along  the  interface  implies  that  the  stresses  are  not  uniform  along  the 


interface  (see  Figure  4).  Along  the  interface,  the  effect  of  dissimil£ir 
material  properties  (i.e.,  G2/G1  <  1)  between  adherend  and  adhesive 
causes  the  normal  stress  to  possess  a  weak  singularity  (Folias  1989) 
near  the  free  surface  z/h=l.  On  the  other  hand,  there  will  be  no  stress 
singularity  present  when  the  material  properties  of  the  adherends  are 
the  same  as  that  of  the  adhesive  (i.e.,  a  uniform  plate). 

In  the  central  region  (see  Figure  4),  the  normal  stress  is 
uniformly  distributed  and  is  slightly  higher  than  the  average  applied 
stress,  and  the  shear  stress  almost  equals  zero  except  in  the  region 
close  to  the  free  surface.  On  the  mid-plane  of  the  adhesive,  the  normal 
stress  decreases  to  a  low  value  near  the  free  surface,  and  the  shear 
stress  is  always  zero  since  this  is  a  symmetrical  problem.  It  is  also 
noted  from  Figure  4  that  there  is  a  significant  variation  in  the  normal 
and  shear  stresses  across  the  adhesive  thickness  especially  near  the 
free  surface  z/h  ==  1  (see  Figure  5).  Near  the  free  surface,  z/h=  0.99, 
the  interface  normal  and  shear  stress  increases  as  the  a/h  value 
increases  (see  Figure  6),  but  it  does  not  change  after  the  a/h  value  is 
greater  than  1.  It  is  also  observed  from  Figure  6  that  the  stress  values 
decreases  as  the  G2/G1  value  Increases. 

In  the  present  analysis,  the  theory  of  linear  elasticity  has  been 
used  to  evaluate  the  stress  field  in  the  adhesive  layer.  In  order  to 
examine  the  possible  failure  mode  of  this  adhesive  joint,  one  may 
choose  the  octahedral  shear  stress,  a  suitable  parameter.  The 

numerical  results  show  that  the  region  around  the  point  of 
intersection  between  the  free  surface  and  the  joint  interface  is 


subjected  to  higher  stress  values.  ThL:s  yielding  is  most  likely  to 
initiate  first  in  this  region.  Figure  7  shows  the  yielding  zone  at  the 
comer  of  the  adhesive  joint  with  respect  for  different  values  of  the 
applied  load.  It  is  observed  that,  when  the  load  is  small  (ie.. 
Oo=0.75aY),  the  initial  yielding  is  along  the  interface.  As  the  load 
increases  further  (ie..  Oo=0.85aY).  this  yielding  area  is  not  extending 
along  the  interface  but  with  an  angle  of  approximately  60°  from  the 
interface.  Thus,  a  crack  will  initiate  at  the  comer  and  propagate  along 
the  interface  up  to  a  certain  distance  beyond  which  it  will  curve  into 
the  adhesive  layer  due  to  the  mixed-mode  (normal  and  shear)  stress 
field  around  the  crack  tip. 

Figure  8  and  9  show  the  comparison  of  the  analytical  results 
obtained  by  this  study  (plane  stress  case)  with  the  results  obtained  by 
the  Finite  Element  Method  (Alwar  et  al.  1976)  with  respect  to  the 
stress  distributions  at  the  interface.  In  this  comparison,  the  a/h  value 
is  0.1.  the  G2/G1  value  is  0.01,  the  Poisson’s  ratio  =  0.3  and  V2  = 

0.33  and  both  results  are  in  fairly  good  agreement.  Except  in  the 
vicinity  of  the  comer  point  z=h  where  a  boundaiy  layer  region  is 
present. 


CONCLUSIONS  AND  RECOMMENDATIONS 

From  the  above  study  of  adhesive  butt  joints,  the  following 


conclusions  can  be  made: 


(1)  The  interface  does  not  remain  plane  under  a  tension  load 
due  to  the  dissimilarity  of  the  material  properties  between  adherend 
and  adhesive,  and  this  displacement  variation  results  in  the  presence 
of  a  stress  singularity  at  the  interface  and  close  to  the  free  surface. 
Along  the  interface,  the  normal  stress  attains  its  maximum  at  the  free 
surface  while  the  shear  stress  attains  its  maximum  very  very  close  to 
the  free  surface.  Furthermore,  it  is  also  observed  that  there  is 
significant  variation  in  the  normal  and  shear  stress  profiles  across  the 
adhesive  thickness  especially  near  the  free  surface. 

(2)  The  maximum  normal  and  shear  stress  increases  as  one 
decreases  the  shear  modulus  of  the  adhesive  layer  or  increases  the  a/h 
value,  but  the  stresses  do  not  change  very  much  when  the  a/h  value  is 
greater  than  1. 

(3)  An  examination  of  the  octahedral  shear  stress  shows  that  a 
crack  is  most  likely  to  initiate  at  the  comer  and  propagate  along  the 
interface  up  to  a  certain  distance  beyond  which  it  will  curve  into  the 
adhesive  layer. 

(4)  The  analytical  results  from  this  study  show  a  good  agreement 
with  the  results  obtained  by  Finite  Element  Method. 

In  this  paper,  the  plane  strain  case  of  adhesive  butt  joint  under 
tension  has  been  solved,  but  the  similar  procedures  can  be  applied  to 
the  problem  of  axisymmetric  butt  Joint  under  tension  or  torsion,  if  one 
transforms  the  general  solution  (16)-(18)  from  rectangular  to 


cylindrical  coordinates.  Furthermore,  the  single-lap  and  bevel  joints 
can  also  be  solved  by  the  application  of  the  general  solution  used  in 
this  analysis.  Identical  and  isotropic  adherends  are  considered  in  this 
present  analysis.  Therefore,  stress  analysis  of  adhesive  joint  with 
nonidentical  or  anisotropic  adherends  is  an  important  extension  of 
this  work.  Finally,  perfect  bonding  is  assumed  in  the  present  analysis, 
but  partial  debonding  (flaw)  is  often  found  at  joint  interface  in 
practical  situations.  Another  important  extension  of  this  work  is  to 
analyze  an  adhesive  joint  which  contains  flaw  at  the  interface  by 
relaxing  the  assumption  of  perfect  bonding. 
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Figure  1  Some  common  engineering  adhesive  joints. 
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Figure  5  Stress  Uxx  and  Xja  for  the  adhesive  layer  across  the 
thickness  at  z/h=  0.98  for  a/h=0.1,  G2/Gi=:0.05, 

Vi=:0.3  and  V2=0.33. 


Vi=0.3.  V2=0.33 


Stress  Oxjc  for  the  adhesive  layer  with  respect  to  the 
different  a/h  values  at  z/h=0.99.  x/ a=  1  for 
VisO.3  and  V2=0.33. 


Figure  6 


r/a  =  0.05 


Figure  7  The  yielding  zone  at  the  comer  of  the  adhesive  joint 
with  respect  to  different  values  of  the  applied  load. 
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and  their  respective  solutions  for  the  displacement  and  stress  fields  are  then  used  in  order  to 
provide  us  with  some  answers  to  the  following  fundamental  questions:  transverse  strength, 
longitudinal  strength,  residual  stresses  due  to  thermal  expansion  mismatch,  modeling  of 
fiber  manix  interface,  edge  effects. 

The  3D  results  are  then  used  to  first  identify  critical  locations  where  failure,  due  to 
fracture,  is  most  likely  to  initiate  and  second  to  derive  fracture  criteria  for  crack  initiation  at 
the  local  level.  The  criteria  reveal  the  dependance  of  the  composite  strength  on  the  material 
properties,  the  local  cell  geometry,  the  ratio  of  the  fiber  volume  fraction,  the  ratio  of  fiber 
radius  to  fiber  length  and  finally  the  applied  mechanical  and  or  thermal  loads. 


